
# The basic libraries
import numpy as np 
import pandas as pd
import os
import math
import scipy

# To plot data
import matplotlib
from matplotlib import pyplot as plt
import matplotlib.image as mpimg
matplotlib.use('Agg')
import sys

# To read in images
import imageio
# To process images
import cv2

# To read ROIs generated in ImageJ
from read_roi import read_roi_file
from read_roi import read_roi_zip

# To fit ellipse
from skimage.measure import EllipseModel
from matplotlib.patches import Ellipse


# Import the required libraries
from numpy import array, savetxt


# Select the txt file (a file dialog will open)
try:
    from tkinter import Tk
    from tkFileDialog import askopenfilenames
except:
    from tkinter import Tk
    from tkinter import filedialog

Tk().withdraw() # we don't want a full GUI, so keep the root window from appearing
filenames = filedialog.askopenfilenames() # show an "Open" dialog box and return the path to the selected file
filename = ''.join(filenames)
print(filename)


# Read in file, and get the number of outlines listed in the file
reader = open(filename,'r')
try:
    #print(reader.read())
    value_lines = list(reader)
finally:
    reader.close()
size_val = len(value_lines)

# Get the number of outlines listed in the file and check, if not correct, adjust it
vals = value_lines[size_val-1].strip("\n")
vals2 = vals.split("\t")
n_outlines = int(vals2[0][-1])+1
print('Number of outlines in this file are:', n_outlines)
#n_outlines = 3
#print(n_outlines)


# Create a new folder to store the separate outlines in
# Get the file path (file_path) and file name (report_name) without the extension
filename_split = filename.split("/")
filename_split_length = len(filename_split)
file_path = (filename[0:(len(filename)-len(filename_split[filename_split_length-1]))])
print(file_path)

report_nametxt = filename_split[filename_split_length-1]
report_name = report_nametxt[0:(len(report_nametxt)-4)]
print(report_name)

# Create a new folder to store the header and the separate outline files in
newpath = '%sFolder_%s' % (file_path,report_name)
print(newpath)
if not os.path.exists(newpath):
    os.makedirs(newpath)


# To save outlines with a file number, starting from a certain number, change count_start to the desired number.
count_start = 10;


# For each outline (j), list the x and y values in separate arrays and write the arrays to a csv file
for j in range(0,n_outlines):
    x_array = np.array([])
    y_array = np.array([])
    
    for i in range(1,size_val):
       # print('i= %i and j= %j',(i,j))
        value_lines_stripped = value_lines[i].strip("\n")
        value_lines_split = value_lines_stripped.split("\t")
    
        if j == (n_outlines-1) and i == (size_val-1): # j is the last outline and i is the last written line 
            # >> append array, save array, save header
            xval = float(value_lines_split[1])
            xvalr = round(xval)
            yval = float(value_lines_split[2])
            yvalr = round(yval)
            x_array = np.append(x_array,xvalr)
            y_array = np.append(y_array,yvalr)
            header_outline = value_lines[0].strip("\n")
            np.savetxt('%s/Outline_%i.csv' % (newpath, j+count_start), np.array([x_array, y_array]).T, fmt='%i', delimiter=',')#, header = header_cont)
            np.savetxt('%s/Header_%s.csv' % (newpath,report_name), [0], header = header_outline)
        elif int(value_lines_split[0][-1]) == j: # if j is not yet 9, and line i belongs to j >> append array
            xval = float(value_lines_split[1])
            xvalr = round(xval)
            yval = float(value_lines_split[2])
            yvalr = round(yval)
            x_array = np.append(x_array,xvalr)
            y_array = np.append(y_array,yvalr)
        elif int(value_lines_split[0][-1]) > j: # if j is not yet 9, and line i is larger than j >> save text
            np.savetxt('%s/Outline_%i.csv' % (newpath, j+count_start), np.array([x_array, y_array]).T, fmt='%i', delimiter=',')#, header = header_cont  
            i = size_val


outline_data = pd.read_csv('Outlines/outlines138', header=None)
n_pixels_outline = round(outline_data.size / 2)
outline_matrix = np.zeros((1024,1024))
for i in range(n_pixels_outline):
    read_x = round(outline_data.loc[i,0])
    read_y = round(outline_data.loc[i,1])
    outline_matrix[read_x,read_y] = 1

# Fill the outline image
filled_binary = scipy.ndimage.binary_fill_holes(outline_matrix)

# Plot shape
%matplotlib inline
plt.imshow(filled_binary, cmap = 'gray')

<matplotlib.image.AxesImage at 0x7f9ea16c3f10>


# Note: some of our images were overlapping each other, which caused some gastruloids to appear twice in our data set.
# We carefully checked all the extracted outlines and removed the doubles. Therefore, 500 - 413 = 87 extracted outlines were removed.

# 96 h gastruloids
a1 = np.arange(1,7)
a2 = np.arange(10,27)
a3 = np.arange(29,37)
a4 = np.array([38, 40])
a5 = np.arange(42,90)
a6 = np.arange(91,100)
a62 = np.array([102,104,106,107,109,112,118,121])
a7 = np.arange(123,157)
a8 = np.array([158,166])
a9 = np.arange(168,172)
a10 = np.array([173])
a11 = np.arange(175,183)
a12 = np.array([185,186,188])
a13 = np.arange(192,214)
a14 = np.arange(216,224)
a15 = np.array([225,226])
a16 = np.arange(228,243)
a17 = np.arange(244,252)
a18 = np.array([253,254])
a19 = np.arange(256,273)
a20 = np.array([274])
a21 = np.arange(276,282)
a22 = np.array([283])
a23 = np.arange(285,296)
a24 = np.array([297])
a25 = np.arange(299,304)
a26 = np.arange(306,311)
a27 = np.array([312,317,319,322,323])
a28 = np.arange(326,349)
a29 = np.arange(350,363)
a30 = np.arange(364,372)
a31 = np.arange(373,376)
a32 = np.arange(377,380)
a33 = np.arange(381,390)
a34 = np.arange(391,400)
a35 = np.arange(401,408)
a36 = np.arange(409,412)
a37 = np.arange(413,435)
a38 = np.arange(437,448)
a39 = np.arange(449,452)
a40 = np.arange(453,459)
a41 = np.arange(461,472)
a42 = np.array([473,474])
a43 = np.arange(476,485)
a44 = np.arange(486,490)
a45 = np.arange(492,500)

a_96h = np.concatenate((a1,a2,a3,a4,a5,a6,a62,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18,a19,a20,a21,a22,a23,a24,a25,a26,a27,a28,a29,a30,a31,a32,a33,a34,a35,a36,a37,a38,a39,a40,a41,a42,a43,a44,a45),dtype = np.int64)
print('Total number of 96 h gastruloids:', len(a_96h))

# 72 h gastruloids
a1 = np.arange(1,11)
a2 = np.arange(13,27)
a3 = np.array([30])
a4 = np.arange(32,89)
a5 = np.arange(91,101)
a6 = np.array([112])
a7 = np.array([114])
a8 = np.arange(116,127)
a9 = np.arange(130,134)
a10 = np.arange(136,140)
a11 = np.arange(146,147)
a12 = np.arange(149,155)
a13 = np.arange(159,171)

a_72h = np.concatenate((a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13))#,dtype = np.int64)
print('Total number of 72 h gastruloids:', len(a_72h))

# 10 uM ROCK inhibitor gastruloids
a1 = np.arange(1,19)
a2 = np.array([20])
a3 = np.arange(23,37)
a4 = np.array([38])
a5 = np.array([40])
a6 = np.arange(42,73)
a7 = np.arange(74,76)
a8 = np.arange(77,112)
a9 = np.array([113])
a10 = np.arange(115,118)
a11 = np.array([119])
a12 = np.array([121])
a13 = np.array([123])

a_10uM_RI = np.concatenate((a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13))#,dtype = np.int64)
print('ROCK inhibitor experiment - total number of treated gastruloids:', len(a_10uM_RI))

# No ROCK inhibitor gastruloids
a1 = np.arange(1,30)
a2 = np.arange(31,34)
a3 = np.arange(35,50)
a4 = np.arange(53,58)
a5 = np.arange(60,67)
a6 = np.arange(68,80)
a7 = np.arange(81,96)
a8 = np.arange(97,102)
a9 = np.arange(103,105)
a10 = np.arange(106,113)
a11 = np.array([115])
a12 = np.arange(117,122)

a_no_RI = np.concatenate((a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12))#,dtype = np.int64)
print('ROCK inhibitor experiment - total number of control gastruloids:', len(a_no_RI))

Total number of 96 h gastruloids: 413
Total number of 72 h gastruloids: 132
ROCK inhibitor experiment - total number of treated gastruloids: 110
ROCK inhibitor experiment - total number of control gastruloids: 106


folder_72h = "/Users/esmeeadegeest/Documents/Outlines_72h_143unique/Outlines"
L2_div_L1_72h_array=[]
L3_div_L1_72h_array=[]
for i in range(np.size(a_72h)):
    i_string = str(a_72h[i])
    data_file = "Contour_"+i_string+"outline.txt"
    path_i = os.path.join(folder_72h,data_file)
    contents = pd.read_table(path_i)
    
    L_values_72h = pd.DataFrame(contents).to_numpy() # convert dataframe to array in numpy
    L2_div_L1_72h = np.divide(L_values_72h[1], L_values_72h[0]) # divide L2 by L1
    L3_div_L1_72h = np.divide(L_values_72h[2], L_values_72h[0]) # divide L3 by L1
    
    L2_div_L1_72h_array.append(L2_div_L1_72h) # store L2/L1 values
    L3_div_L1_72h_array.append(L3_div_L1_72h) # store L3/L1 values


folder_96h = "/Users/esmeeadegeest/Documents/L_values_no_doubles_413/"
L2_div_L1_96h_array=[]
L3_div_L1_96h_array=[]
for i in range(np.size(a_96h)):
    i_string = str(a_96h[i])
    data_file = "Data_"+i_string+".txt"
    path_i = os.path.join(folder_96h,data_file)
    contents = pd.read_table(path_i)
    
    L_values_96h = pd.DataFrame(contents).to_numpy() # convert dataframe to array in numpy
    L2_div_L1_96h = np.divide(L_values_96h[1], L_values_96h[0]) # divide L2 by L1
    L3_div_L1_96h = np.divide(L_values_96h[2], L_values_96h[0]) # divide L3 by L1
    
    L2_div_L1_96h_array.append(L2_div_L1_96h) # store L2/L1 values
    L3_div_L1_96h_array.append(L3_div_L1_96h) # store L3/L1 values


folder_No_RI = "/Users/esmeeadegeest/Documents/L-values_No-RI/"
L2_div_L1_No_RI_array=[]
L3_div_L1_No_RI_array=[]
for i in range(np.size(a_no_RI)):
    i_string = str(a_no_RI[i])
    data_file = "Contour_"+i_string+"outline.txt"
    path_i = os.path.join(folder_No_RI,data_file)
    contents = pd.read_table(path_i)
    
    L_values_No_RI = pd.DataFrame(contents).to_numpy() # convert dataframe to array in numpy
    L2_div_L1_No_RI = np.divide(L_values_No_RI[1], L_values_No_RI[0]) # divide L2 by L1
    L3_div_L1_No_RI = np.divide(L_values_No_RI[2], L_values_No_RI[0]) # divide L3 by L1
    
    L2_div_L1_No_RI_array.append(L2_div_L1_No_RI) # store L2/L1 values
    L3_div_L1_No_RI_array.append(L3_div_L1_No_RI) # store L3/L1 values


folder_10uM_RI = "/Users/esmeeadegeest/Documents/L-values_10uM-RI/"
L2_div_L1_10uM_RI_array=[]
L3_div_L1_10uM_RI_array=[]
for i in range(np.size(a_10uM_RI)):
    i_string = str(a_10uM_RI[i])
    data_file = "Contour_"+i_string+"outline.txt"
    path_i = os.path.join(folder_10uM_RI,data_file)
    contents = pd.read_table(path_i)
    
    L_values_10uM_RI = pd.DataFrame(contents).to_numpy() # convert dataframe to array in numpy
    L2_div_L1_10uM_RI = np.divide(L_values_10uM_RI[1], L_values_10uM_RI[0]) # divide L2 by L1
    L3_div_L1_10uM_RI = np.divide(L_values_10uM_RI[2], L_values_10uM_RI[0]) # divide L3 by L1
    
    L2_div_L1_10uM_RI_array.append(L2_div_L1_10uM_RI) # store L2/L1 values
    L3_div_L1_10uM_RI_array.append(L3_div_L1_10uM_RI) # store L3/L1 values


#72 h and 96 h
%matplotlib inline
font = {'size'   : 12}
matplotlib.rc('font', **font)
matplotlib.rcParams['font.sans-serif'] = "Arial"
matplotlib.rcParams['font.family'] = "sans-serif" 
cbmap=['#377eb8', '#ff7f00', '#4daf4a',
                  '#f781bf', '#a65628', '#984ea3',
                  '#999999', '#e41a1c', '#dede00']

plt.gcf().clear()
fig = plt.figure(1)
ax = fig.add_subplot(111)

ax.scatter(L2_div_L1_72h_array,L3_div_L1_72h_array, s=1, marker ="o", label='72h gastruloids', c=cbmap[5])
ax.scatter(L2_div_L1_96h_array,L3_div_L1_96h_array, s=1, marker ="o", label='96h gastruloids', c=cbmap[1])

fig.set_size_inches(10/2.54, 8/2.54)
plt.xlabel('$L_2/L_1$')
plt.ylabel('$L_3/L_1$')
#plt.title('all <-> all')
handles, labels = ax.get_legend_handles_labels()
lgd = ax.legend(handles, labels, loc='upper left', bbox_to_anchor=(1,1))
fig.set_size_inches(10/2.54, 8/2.54)

# save figure
#fig.savefig('Figures_paper/Figure_1_LOCO-EFA_72h_96h_purple_orange.svg', bbox_extra_artists=([lgd]), bbox_inches='tight')


# ROCK inhibitor experiment
%matplotlib inline
font = {'size'   : 12}
matplotlib.rc('font', **font)
matplotlib.rcParams['font.sans-serif'] = "Arial"
matplotlib.rcParams['font.family'] = "sans-serif" 
cbmap=['#377eb8', '#ff7f00', '#4daf4a',
                  '#f781bf', '#a65628', '#984ea3',
                  '#999999', '#e41a1c', '#dede00']

plt.gcf().clear()
fig = plt.figure(1)
ax = fig.add_subplot(111)

ax.scatter(L2_div_L1_No_RI_array,L3_div_L1_No_RI_array, s=1, marker ="o", label='Control', c=cbmap[1])
ax.scatter(L2_div_L1_10uM_RI_array,L3_div_L1_10uM_RI_array, s=1, marker ="o", label=r'10 $\mu$M Y-27632', c=cbmap[5])

fig.set_size_inches(10/2.54, 8/2.54)
plt.xlabel('$L_2/L_1$')
plt.ylabel('$L_3/L_1$')
plt.xlim(-0.02, 0.7)
plt.ylim(-0.015, 0.305)
#plt.title('all <-> all')
handles, labels = ax.get_legend_handles_labels()
lgd = ax.legend(handles, labels, loc='upper left', bbox_to_anchor=(1,1))
fig.set_size_inches(10/2.54, 8/2.54)

# save figure
#fig.savefig('Documents/Figure_X_LOCO-EFA_RI_purple_orange.svg', bbox_extra_artists=([lgd]), bbox_inches='tight')


# Pre-processing in FIJI:
# run("RGB Color");
# run("Enhance Contrast", "saturated=0.35");
# should be stored as PNG file.)

%matplotlib inline
img = []
img_GB = []
img_GB_thresh = []
img_binary = []

plt.rcParams['figure.figsize'] = [30, 20]
kernel = np.ones((5, 5), np.uint8) # Taking a matrix of size 5 as the kernel for the morphological operations

# Perform morphological operations on the maximum z projection for each time point using a specified threshold value and store the images in img_binary
for i in range(1,32):
    img_read = cv2.imread('Mosaic_TL_mCherry/Max_projections_Z/MAX_T'+str(i)+'_TL-P1-10min-33ts_G2_561nm-1700mA_160ms_6p7_4um_40zimage1_crop.png',0)
    img.append(img_read)
    
    if i == 1:
        threshold_value = 63 # We used 40 for all frames, except for frame of t = 1 (63) and t = 26 (45).
    elif i == 26:
        threshold_value = 45
    else:
        threshold_value = 40
        
    # Apply Gaussian blur
    img_read_GB = cv2.GaussianBlur(img_read, (51,51), cv2.BORDER_DEFAULT) # kernel size should be odd
    img_GB.append(img_read_GB)
    # Apply threshold value >> check in the next cell for each time frame whether the threshold value is correct
    ret,img_GB_read_thresh1 = cv2.threshold(img_read_GB,threshold_value,255,cv2.THRESH_BINARY)
    img_GB_thresh.append(img_GB_read_thresh1)
    # Apply dilation multiple times to fill the holes. 
    img_dilation = cv2.dilate(img_GB_read_thresh1, kernel, iterations=40)
    # Apply subsequently erosion to retrieve the initial size
    img_erosion = cv2.erode(img_dilation, kernel, iterations=40)
    img_binary.append(img_erosion)


%matplotlib inline
plt.rcParams['figure.figsize'] = [15, 15]
fig, axs = plt.subplots(2, 2)

# Show for one time point the original image, the image after Gaussian blur, 
# the image after applying a threshold, and the image after dilation and erosion.
n = 31 # time point
axs[0,0].imshow(img[n-1], cmap='gray')
axs[0,1].imshow(img_GB[n-1], cmap='gray')
axs[1,0].imshow(img_GB_thresh[n-1], cmap='gray')
axs[1,1].imshow(img_binary[n-1], cmap='gray')

plt.show()


# Save binary images, note down the threshold value for each image
# We used threshold value 40 for all frames, except for the frame of t = 1 (63) and t = 26 (45).

# t = 1, threshold = 63
cv2.imwrite('Mosaic_TL_mCherry/Binary_images/t1-63-binary.png', img_binary[0])

# t = 2-25, threshold = 40
for i in range(2,26):
    cv2.imwrite('Mosaic_TL_mCherry/Binary_images/t'+str(i)+'-40-binary.png', img_binary[i-1])
    
# t = 26, threshold = 45
cv2.imwrite('Mosaic_TL_mCherry/Binary_images/t26-45-binary.png', img_binary[25])

# t = 27-31, threshold = 40
for i in range(27,32):
    cv2.imwrite('Mosaic_TL_mCherry/Binary_images/t'+str(i)+'-40-binary.png', img_binary[i-1])


# Area calculation and comparison t1 and t31
image_loc_t1 = 'Mosaic_TL_mCherry/Binary_images/t1-63-binary.png'
image_loc_t31 = 'Mosaic_TL_mCherry/Binary_images/t31-40-binary.png'

mask_t1 = cv2.imread(image_loc_t1)
mask_t31 = cv2.imread(image_loc_t31)

imgray_t1 = cv2.cvtColor(mask_t1,cv2.COLOR_BGR2GRAY)
imgray_t31 = cv2.cvtColor(mask_t31,cv2.COLOR_BGR2GRAY)

area_t1 = np.sum(imgray_t1)/255
area_t31 = np.sum(imgray_t31)/255
print('Ratio area t31/t1:', area_t31/area_t1)

plt.rcParams['figure.figsize'] = [10, 10]
fig, axs = plt.subplots(1, 2)
axs[0].imshow(mask_t1, cmap='gray')
axs[1].imshow(mask_t31, cmap='gray')

Ratio area t31/t1: 1.3836912100209595

<matplotlib.image.AxesImage at 0x7fb20addab80>


# We used the following code to extract outlines from the binary images.
def extract_outline_iv():

    for j in range(1,32):
        ind = list() # index
        cont_all = list() # contours
        
        if j == 1:
            image_loc = 'Mosaic_TL_mCherry/Binary_images/t'+str(j)+'-63-binary.png'
            file_loc = 'Mosaic_TL_mCherry/Binary_outlines/t'+str(j)+'-63-binary_outline.txt'
        elif j == 26:
            image_loc = 'Mosaic_TL_mCherry/Binary_images/t'+str(j)+'-45-binary.png'
            file_loc = 'Mosaic_TL_mCherry/Binary_outlines/t'+str(j)+'-45-binary_outline.txt'
        else:
            image_loc = 'Mosaic_TL_mCherry/Binary_images/t'+str(j)+'-40-binary.png'
            file_loc = 'Mosaic_TL_mCherry/Binary_outlines/t'+str(j)+'-40-binary_outline.txt'
            
        file1 = open(file_loc,"a")
        mask = cv2.imread(image_loc)

        if mask is not None:

            imgray = cv2.cvtColor(mask,cv2.COLOR_BGR2GRAY)
            ret,thresh = cv2.threshold(imgray,200,255,0)
            contours, im2 = cv2.findContours(thresh,cv2.RETR_TREE,cv2.CHAIN_APPROX_NONE)         
            maxsizeloc = 0
            maxsize = 0
           
            for i in range (len(contours)):
                if contours[i].shape[0] > maxsize:
                    maxsizeloc = i
                    maxsize = contours[i].shape[0]
            cont = np.zeros((contours[maxsizeloc].shape[0], 2))
            cont[:,0] = contours[maxsizeloc][:,0,0]
            cont[:,1] = contours[maxsizeloc][:,0,1]

            plt.plot(cont[:,0], cont[:,1])
            cont_all.append(cont)
            
            # Write outlines in txt file
            for i in range (maxsize):
                file1.write(str(cont_all[0][i][0]) + ',' + str(cont_all[0][i][1]) + '\n')
            file1.close()
          
    return(cont_all)

extract_outline_iv()


# Ellipse fit parameters
xc_iv = np.empty(31) # center coordinate x of ellipse
yc_iv = np.empty(31) # center coordinate y of ellipse
a_iv = np.empty(31) # length short arm of ellipse
b_iv = np.empty(31) # length long arm of ellipse
theta_iv= np.empty(31) # rotation angle of ellipse

# Figure
%matplotlib inline
fig, axs = plt.subplots(2, 16,figsize=(100,30))#sharex=True, sharey=True)

# Load outline coordinates for each time step and fit ellipse, store ellipse parameters
for j in range (1,32):
    # txt file name in which outline is stored 
    if j == 1:
        outline_loc = 'Mosaic_TL_mCherry/Binary_outlines/t'+str(j)+'-63-binary_outline.txt'
        file_loc = 'Mosaic_TL_mCherry/Binary_outlines/t'+str(j)+'-63_xycenter-theta-axes.txt' # assign a file location and open it to later on write the ellipse parameters in the file
    elif j == 26:
        outline_loc = 'Mosaic_TL_mCherry/Binary_outlines/t'+str(j)+'-45-binary_outline.txt'
        file_loc = 'Mosaic_TL_mCherry/Binary_outlines/t'+str(j)+'-45_xycenter-theta-axes.txt'
    else:
        outline_loc = 'Mosaic_TL_mCherry/Binary_outlines/t'+str(j)+'-40-binary_outline.txt'
        file_loc = 'Mosaic_TL_mCherry/Binary_outlines/t'+str(j)+'-40_xycenter-theta-axes.txt'
   
    file1 = open(file_loc,"a")
    
    # open outline file and store coordinates in array
    with open(outline_loc) as file_name:
        points = np.loadtxt(file_name, delimiter=",")
    a_points = np.array(points)
    x = a_points[:, 0]
    y = a_points[:, 1]
    
    # fit ellipse to array and get parameters
    ell = EllipseModel()
    ell.estimate(a_points)
    xc_iv[j-1], yc_iv[j-1], a_iv[j-1], b_iv[j-1], theta_iv[j-1] = ell.params
    
    # Write parameters to txt file
    #file1.write(str(xc_iv[j-1]) + ',' + str(yc_iv[j-1]) + ',' + str(theta_iv[j-1]) + ',' + str(a_iv[j-1]) + ',' + str(b_iv[j-1]) + '\n')
    
    # Plot ellipse fits
    if j < 16: # top row of images
        axs[0,j-1].scatter(x, y) # outline
        axs[0,j-1].scatter(xc_iv[j-1], yc_iv[j-1], color='red', s=100) # center of ellipse
        axs[0,j-1].set_xlim(0,1600)
        axs[0,j-1].set_ylim(0,3000)
        axs[0,j-1].set_xticks([])
        axs[0,j-1].set_yticks([])
        
        # ellipse
        ell_patch = Ellipse((xc_iv[j-1], yc_iv[j-1]), 2*a_iv[j-1], 2*b_iv[j-1], theta_iv[j-1]*180/np.pi, edgecolor='red', facecolor='none')
        # ellipse with orientation angle 0
        ell_patch2 = Ellipse((xc_iv[j-1], yc_iv[j-1]), 2*a_iv[j-1], 2*b_iv[j-1], 0, edgecolor='black', facecolor='none')

        axs[0,j-1].add_patch(ell_patch)
        axs[0,j-1].add_patch(ell_patch2)
        
    else: # bottom row of images
        axs[1,j-16].scatter(x, y) # outline
        axs[1,j-16].scatter(xc_iv[j-1], yc_iv[j-1], color='red', s=100) # center of ellipse
        axs[1,j-16].set_xlim(0,1600)
        axs[1,j-16].set_ylim(0,3000)
        axs[1,j-16].set_xticks([])
        axs[1,j-16].set_yticks([])
        
        # ellipse
        ell_patch = Ellipse((xc_iv[j-1], yc_iv[j-1]), 2*a_iv[j-1], 2*b_iv[j-1], theta_iv[j-1]*180/np.pi, edgecolor='red', facecolor='none')
        # ellipse with orientation angle 0
        ell_patch2 = Ellipse((xc_iv[j-1], yc_iv[j-1]), 2*a_iv[j-1], 2*b_iv[j-1], 0, edgecolor='black', facecolor='none')

        axs[1,j-16].add_patch(ell_patch)
        axs[1,j-16].add_patch(ell_patch2)
plt.show()

# From the images we can tell that for some shapes the short arm length and long arm length are
# stored in the wrong parameters: b and a, respectively, while they should be stored in a and b, respectively.


# Here we print the short arm length (a) and long arm length (b) and see that they are not always stored correctly.
print('Time point 1 (correct): a = ',a_iv[0],'b =', b_iv[0])
print('Time point 4 (incorrect): a = ',a_iv[3],'b = ', b_iv[3])

Time point 1 (correct): a =  588.7241640888858 b = 879.3560145828596
Time point 4 (incorrect): a =  921.3089122428711 b =  591.6155040634976


# We traced cells manually using the point tool in ImageJ. The trajectories were stored in an ROI set.
# Import the roi_set of each cell
N_cells = 14
cell_list = [None] * N_cells
roi_set = [None] * N_cells

for i in range(N_cells):
    path_file = 'Mosaic_TL_mCherry/Cells_by_number/RoiSet_Cell_' + str(i+1) + '.zip'
    roi_set[i] = read_roi_zip(path_file)


# Dividing cell script for cell 3 (alternating position)
n = 3
roi_set_sel = roi_set[n-1]
time_step_division = 17
t = 1 # start time

division_status = 'Dividing cell'
path_x = np.arange(0) # mother cell x
path_y = np.arange(0) # mother cell y
path_x_d1 = np.arange(0) # daughter cell 1 x
path_x_d2 = np.arange(0) # daughter cell 2 x
path_y_d1 = np.arange(0) # daughter cell 1 y
path_y_d2 = np.arange(0) # daughter cell 2 xy
roi_set_length = len(roi_set_sel)
d1 = bool(True) # daughter cell 1 is set to true

for key, value in roi_set_sel.items():
    path_roi = roi_set_sel[key]
    
    if t < time_step_division:
        path_x = np.append(path_x, path_roi['x'][0])
        path_y = np.append(path_y, path_roi['y'][0])
        t = t + 1
    elif t >= time_step_division:
        if d1 == True:
            path_x_d1 = np.append(path_x_d1, path_roi['x'][0])
            path_y_d1 = np.append(path_y_d1, path_roi['y'][0])
            d1 = bool(False) # set daughter cell 1 to false > store next position in array of 2nd daughter cell
            t = t + 1
        elif d1 == False:
            path_x_d2 = np.append(path_x_d2, path_roi['x'][0])
            path_y_d2 = np.append(path_y_d2, path_roi['y'][0])
            d1 = bool(True) # set daughter cell 1 to true > store next position in array of 1st daughter cell
            t = t + 1
            
# store the trajectories in a list   
cell_list[n-1] = [path_x, path_y, division_status, time_step_division, path_x_d1[:(len(path_x_d1))], path_y_d1[:(len(path_y_d1))], path_x_d2[:(len(path_x_d2))], path_y_d2[:(len(path_y_d2))]]


# Dividing cell script for cell 5 (weird order: first alternating positions, then paths subsequently stored in ROI set)
n = 5
roi_set_sel = roi_set[n-1]
time_step_division_par = 9
time_step_division_sep = 16
t_sep_start = time_step_division_par + 2*(time_step_division_sep - time_step_division_par)
t_sep_d2_start = t_sep_start + (31 - time_step_division_sep + 1)

roi_set_length = len(roi_set_sel)
t = 1
division_status = 'Dividing cell'
path_x = np.arange(0)
path_y = np.arange(0)
path_x_d1 = np.arange(0)
path_x_d2 = np.arange(0)
path_y_d1 = np.arange(0)
path_y_d2 = np.arange(0)
d1 = bool(True)

for key, value in roi_set_sel.items():
    path_roi = roi_set_sel[key]
    
    if t < time_step_division_par: # to 8 (incl)
        path_x = np.append(path_x, path_roi['x'][0])
        path_y = np.append(path_y, path_roi['y'][0])
        t = t + 1

    elif t >= time_step_division_par and t < t_sep_start: # 9 to 15 (incl)
        if d1 == True:
            path_x_d1 = np.append(path_x_d1, path_roi['x'][0])
            path_y_d1 = np.append(path_y_d1, path_roi['y'][0])
            d1 = bool(False)
            t = t + 1
        elif d1 == False:
            path_x_d2 = np.append(path_x_d2, path_roi['x'][0])
            path_y_d2 = np.append(path_y_d2, path_roi['y'][0])
            d1 = bool(True)
            t = t + 1
    elif t >= t_sep_start and t < t_sep_d2_start:
        path_x_d1 = np.append(path_x_d1, path_roi['x'][0])
        path_y_d1 = np.append(path_y_d1, path_roi['y'][0])
        t = t + 1
    elif t >= t_sep_d2_start:
        path_x_d2 = np.append(path_x_d2, path_roi['x'][0])
        path_y_d2 = np.append(path_y_d2, path_roi['y'][0])
        t = t + 1

cell_list[n-1] = [path_x, path_y, division_status, time_step_division_par, path_x_d1[:(len(path_x_d1))], path_y_d1[:(len(path_y_d1))], path_x_d2[:(len(path_x_d2))], path_y_d2[:(len(path_y_d2))]]


# Dividing cell script for cell 7
n = 7
time_step_division = 14
roi_set_sel = roi_set[n-1]
roi_set_length = len(roi_set_sel)

t = 1
division_status = 'Dividing cell'
path_x = np.arange(0)
path_y = np.arange(0)
path_x_d1 = np.arange(0)
path_x_d2 = np.arange(0)
path_y_d1 = np.arange(0)
path_y_d2 = np.arange(0)

for key, value in roi_set_sel.items():
    path_roi = roi_set_sel[key]
    if t < time_step_division:
        path_x = np.append(path_x, path_roi['x'][0])
        path_y = np.append(path_y, path_roi['y'][0])
        t = t + 1
    elif t >= time_step_division and t <= 31:
        path_x_d1 = np.append(path_x_d1, path_roi['x'][0])
        path_y_d1 = np.append(path_y_d1, path_roi['y'][0])
        t = t + 1
    elif t > 31:
        path_x_d2 = np.append(path_x_d2, path_roi['x'][0])
        path_y_d2 = np.append(path_y_d2, path_roi['y'][0])
        t = t + 1

cell_list[n-1] = [path_x, path_y, division_status, time_step_division, path_x_d1[:(len(path_x_d1))], path_y_d1[:(len(path_y_d1))], path_x_d2[:(len(path_x_d2))], path_y_d2[:(len(path_y_d2))]]


# Dividing cell script for cell 9
n = 9
time_step_division = 12
roi_set_sel = roi_set[n-1]

roi_set_length = len(roi_set_sel)
t = 1
division_status = 'Dividing cell'
path_x = np.arange(0)
path_y = np.arange(0)
path_x_d1 = np.arange(0)
path_x_d2 = np.arange(0)
path_y_d1 = np.arange(0)
path_y_d2 = np.arange(0)

for key, value in roi_set_sel.items():
    path_roi = roi_set_sel[key]
    
    if t < time_step_division:
        path_x = np.append(path_x, path_roi['x'][0])
        path_y = np.append(path_y, path_roi['y'][0])
        t = t + 1

    elif t >= time_step_division and t <= 31:
        path_x_d1 = np.append(path_x_d1, path_roi['x'][0])
        path_y_d1 = np.append(path_y_d1, path_roi['y'][0])
        t = t + 1

    elif t > 31:
        path_x_d2 = np.append(path_x_d2, path_roi['x'][0])
        path_y_d2 = np.append(path_y_d2, path_roi['y'][0])
        t = t + 1

cell_list[n-1] = [path_x, path_y, division_status, time_step_division, path_x_d1[:(len(path_x_d1))], path_y_d1[:(len(path_y_d1))], path_x_d2[:(len(path_x_d2))], path_y_d2[:(len(path_y_d2))]]


# Non-dividing cells: 1, 2, 4, 6, 10, 11, 12, 13, 14
# Fill in here:
n = [1, 2, 4, 6, 10, 11, 12, 13, 14]

division_status = 'Non-dividing cell'

for i in n:
    path_x = np.arange(0)
    path_y = np.arange(0)
    roi_set_sel = roi_set[i-1]
    
    for key, value in roi_set_sel.items():
        path_roi = roi_set_sel[key]
        path_x = np.append(path_x,path_roi['x'][0])
        path_y = np.append(path_y,path_roi['y'][0])
        
    cell_list[i-1] = [path_x[:(len(path_x))], path_y[:(len(path_y))], division_status]


# Non-dividing cell - cell 8 (reverse order)
# Fill in here:
n = 8
roi_set_sel = roi_set[n-1]

division_status = 'Non-dividing cell'
path_x = np.arange(0)
path_y = np.arange(0)

for key, value in roi_set_sel.items():
    path_roi = roi_set_sel[key]
    path_x = np.append(path_x,path_roi['x'][0])
    path_y = np.append(path_y,path_roi['y'][0])

path_x_rev = path_x[::-1]
path_y_rev = path_y[::-1]

cell_list[n-1] = [path_x_rev[:(len(path_x_rev))], path_y_rev[:(len(path_y_rev))], division_status]


n_cells = 14 
divcells = 4 # 4 dividing cells
#cells_sel = [1, 2, 4, 6, 8, 10, 11, 12, 13, 14] #non-dividing cells only
cells_sel = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14] #all cells
div_cells_no = [3, 5, 7, 9]

t_steps = 31
go_twice = False # is true for dividing cells
count = 0
theta_corr_iv = [None] * t_steps

rotated_cell_centers_iv = np.empty((n_cells+divcells, t_steps, 4)) # we store: x, y, V_x, V_y

for i in range(n_cells):

    sel_cell = cells_sel[i]-1
        
    if cell_list[sel_cell][2] == 'Dividing cell':
       
        go_twice = True
    
        x_d1 = []
        y_d1 = []
        x_d2 = []
        y_d2 = []

        x_d1 = np.append([cell_list[sel_cell][0]],[cell_list[sel_cell][4]]) # paste trajectory (x coordinate) daughter cell 1 after trajectory mother cell
        y_d1 = np.append([cell_list[sel_cell][1]],[cell_list[sel_cell][5]]) # paste trajectory (y coordinate) daughter cell 1 after trajectory mother cell
    
        x_d2 = np.append([cell_list[sel_cell][0]],[cell_list[sel_cell][6]]) # paste trajectory (x coordinate) daughter cell 2 to trajectory mother cell
        y_d2 = np.append([cell_list[sel_cell][1]],[cell_list[sel_cell][7]]) # paste trajectory (y coordinate) daughter cell 2 to trajectory mother cell
    
    else:
        go_twice = False
    
    if go_twice == False: # non-dividing cell

        for j in range(t_steps):

            if theta_iv[j] > 1.5708: # if gastruloid (ellipse fit) orientation angle is larger than pi/2, correct the angle by subtracting pi/2
                theta_corr_iv[j] = theta_iv[j] - 1.5708
            else:
                theta_corr_iv[j] = theta_iv[j]
            
            # 1. Calculate the difference in x and y position of a cell w.r.t. the center of the ellipse for each time step
            dx = cell_list[sel_cell][0][j]-xc_iv[j]
            dy = cell_list[sel_cell][1][j]-yc_iv[j]
            
            # 2. Calculate the difference in x and y position of each cell w.r.t. to the new x and y axes, the axes of the ellipse, where the long axis is 
            # the new y axis and the short axis is the new x axis. We call them V_x and V_y.
            V_x = np.cos(theta_corr_iv[j]) * dx + np.sin(theta_corr_iv[j]) * dy
            V_y = np.cos(theta_corr_iv[j]) * dy - np.sin(theta_corr_iv[j]) * dx
            
            # 3. Store V_x and V_y and store the new x and y coordinates: V_x + x(center of ellipse) and V_y + y(center of ellipse)
            rotated_cell_centers_iv[i+count,j,0] = V_x + xc_iv[j]
            rotated_cell_centers_iv[i+count,j,2] = V_x
            
            rotated_cell_centers_iv[i+count,j,1] = V_y + yc_iv[j]
            rotated_cell_centers_iv[i+count,j,3] = V_y
            
                
    elif go_twice == True: # dividing cell
        
        for j in range(t_steps):
            
            if theta_iv[j] > 1.5708: # if gastruloid (ellipse fit) orientation angle is larger than pi/2, correct the angle by subtracting pi/2
                theta_corr_iv[j] = theta_iv[j] - 1.5708
            else:
                theta_corr_iv[j] = theta_iv[j]
              
            # 1. Calculate the difference in x and y position of a cell w.r.t. the center of the ellipse for each time step
            dx = x_d1[j]-xc_iv[j]
            dy = y_d1[j]-yc_iv[j]
            
            # 2. Calculate the difference in x and y position of each cell w.r.t. to the new x and y axes, the axes of the ellipse, where the long axis is 
            # the new y axis and the short axis is the new x axis. We call them V_x and V_y.
            V_x = np.cos(theta_corr_iv[j]) * dx + np.sin(theta_corr_iv[j]) * dy
            V_y = np.cos(theta_corr_iv[j]) * dy - np.sin(theta_corr_iv[j]) * dx
            
            # 3. Store V_x and V_y and store the new x and y coordinates: V_x + x(center of ellipse) and V_y + y(center of ellipse)
            rotated_cell_centers_iv[i+count,j,0] = V_x + xc_iv[j]
            rotated_cell_centers_iv[i+count,j,2] = V_x
            
            rotated_cell_centers_iv[i+count,j,1] = V_y + yc_iv[j]
            rotated_cell_centers_iv[i+count,j,3] = V_y
        
        count = count + 1 # for 2nd daughter cell
        
        for j in range(t_steps):
            
            if theta_iv[j] > 1.5708: # if gastruloid (ellipse fit) orientation angle is larger than pi/2, correct the angle by subtracting pi/2
                theta_corr_iv[j] = theta_iv[j] - 1.5708
            else:
                theta_corr_iv[j] = theta_iv[j]
            
            # 1. Calculate the difference in x and y position of a cell w.r.t. the center of the ellipse for each time step
            dx = x_d2[j]-xc_iv[j]
            dy = y_d2[j]-yc_iv[j]
            
            # 2. Calculate the difference in x and y position of each cell w.r.t. to the new x and y axes, the axes of the ellipse, where the long axis is 
            # the new y axis and the short axis is the new x axis. We call them V_x and V_y.
            V_x = np.cos(theta_corr_iv[j]) * dx + np.sin(theta_corr_iv[j]) * dy
            V_y = np.cos(theta_corr_iv[j]) * dy - np.sin(theta_corr_iv[j]) * dx
            
            # 3. Store V_x and V_y and store the new x and y coordinates: V_x + x(center of ellipse) and V_y + y(center of ellipse)
            rotated_cell_centers_iv[i+count,j,0] = V_x + xc_iv[j]
            rotated_cell_centers_iv[i+count,j,2] = V_x
            
            rotated_cell_centers_iv[i+count,j,1] = V_y + yc_iv[j]
            rotated_cell_centers_iv[i+count,j,3] = V_y


# T = 1, 14 cells starting points

%matplotlib inline
plt.rcParams['figure.figsize'] = [30, 20]
imgplot= mpimg.imread('Mosaic_TL_mCherry/MAX_prediction-t1-200epochsRGB_SB100um.png')
plt.imshow(imgplot)#, extent=[0, 300, 0, 300])

# Colors
cmap = [#'#000000', #black
        '#e69d00', #orange
        '#f0e442', #yellow
        '#009e73', #bluish green
        '#56b4e9', #sky blue
        '#cc79a7', #reddish purple
        '#d55e00', #vermillion (reddish)
        '#0072b2', #blue
        #'#ffffff', #white
        '#7f7f7f'] #grey
cell_marker = ['$1$','$2$','$3$','$4$','$5$','$6$','$7$','$8$','$9$','$10$','$11$','$12$','$13$','$14$']
n = 14

for i in range(n):
    if i >= 9:
        multiplying_factor = 1.5 # use a larger icon for cells that start to count from 10
    else:
        multiplying_factor = 1
       
    i_sel = i
    cell_marker_i = cell_marker[i_sel]
    if i < 8:
        color_i = cmap[i]
    if i >= 8 and i < 16:
        color_i = cmap[i-8] 

    # Plot starting points 
    plt.plot(cell_list[i_sel][0][0],cell_list[i_sel][1][0], color = 'k', marker='o',
     linewidth=2, markersize=23) # start point black edge
    plt.plot(cell_list[i_sel][0][0],cell_list[i_sel][1][0], color = color_i, marker='o',
     linewidth=2, markersize=21) # start point colored
    plt.plot(cell_list[i_sel][0][0],cell_list[i_sel][1][0], color = 'w', marker='o',
     linewidth=2, markersize=16) # start point white
    plt.plot(cell_list[i_sel][0][0],cell_list[i_sel][1][0], color = 'k', marker=cell_marker_i,
     linewidth=2, markersize=10*multiplying_factor) # start point number
          
#get current axes
ax = plt.gca()

#hide x-axis
ax.get_xaxis().set_visible(False)

#hide y-axis 
ax.get_yaxis().set_visible(False)

# Save figure
#plt.savefig("Mosaic_TL_mCherry/Figure_in_vitro_tracking_mosaic_icons_t1_more_white.png", dpi=600, bbox_inches='tight',pad_inches = 0)

# Plot figure
fig = plt.show()


# T = 31, 14 cells trajectories
%matplotlib inline
plt.rcParams['figure.figsize'] = [30, 20]

imgplot= mpimg.imread('Mosaic_TL_mCherry/MAX_prediction-t31-200epochsRGB-SB100um.png')
plt.imshow(imgplot)#, extent=[0, 300, 0, 300])

# Colors
cmap = [#'#000000', #black
        '#e69d00', #orange
        '#f0e442', #yellow
        '#009e73', #bluish green
        '#56b4e9', #sky blue
        '#cc79a7', #reddish purple
        '#d55e00', #vermillion (reddish)
        '#0072b2', #blue
        #'#ffffff', #white
        '#7f7f7f'] #grey
cell_marker = ['$1$','$2$','$3$','$4$','$5$','$6$','$7$','$8$','$9$','$10$','$11$','$12$','$13$','$14$']
n = 14

for i in range(n):
    if i >= 9:
        multiplying_factor = 1.5 # use a larger icon for cells that start to count from 10
    else:
        multiplying_factor = 1
       
    i_sel = i
    cell_marker_i = cell_marker[i_sel]
    if i < 8:
        color_i = cmap[i]
    if i >= 8 and i < 16:
        color_i = cmap[i-8]
    
    # dividing cells
    # division points = cell 3: 17, cell 5: 16, cell 7: 14, cell 9: 12
    if cell_list[i_sel][2] == 'Dividing cell':
        cell_division_point = cell_list[i_sel][3]
        
        plt.plot(cell_list[i_sel][4], cell_list[i_sel][5], color = 'k', #linestyle = 'dotted',
         linewidth = 3, markersize = 3) # black trajectory "daughter cell 1"
        plt.plot(cell_list[i_sel][4], cell_list[i_sel][5], color = color_i,
          linewidth=2.5, markersize=2.5) # colored trajectory "daughter cell 1"
        
        plt.plot(cell_list[i_sel][6], cell_list[i_sel][7], color = 'k', #linestyle = 'dotted',
         linewidth = 3, markersize = 3) # black trajectory "daughter cell 2"
        plt.plot(cell_list[i_sel][6], cell_list[i_sel][7], color = color_i,
          linewidth=2.5, markersize=2.5) # colored trajectory "daughter cell 2"
    
    # Trajectory from start to end (mother cell or non-dividing cell)
    plt.plot(cell_list[i_sel][0], cell_list[i_sel][1], color = 'k', #linestyle = 'dotted',
     linewidth = 3, markersize = 3) # line connecting start and end point black 
    plt.plot(cell_list[i_sel][0], cell_list[i_sel][1], color = color_i,
     linewidth=2.5, markersize=2.5) # line connecting start and end point color

    # start points
    if cell_list[i_sel][2] == 'Dividing cell':
        
        # x and y coordinates for connecting line mother cell and daughter cell
        mother_x = cell_list[i_sel][0][-1]
        mother_y = cell_list[i_sel][1][-1]
        daughter_1_x = cell_list[i_sel][4][0]
        daughter_1_y = cell_list[i_sel][5][0]
        daughter_2_x = cell_list[i_sel][6][0]
        daughter_2_y = cell_list[i_sel][7][0]
        x_points_d1 = [mother_x, daughter_1_x]
        x_points_d2 = [mother_x, daughter_2_x]
        y_points_d1 = [mother_y, daughter_1_y]
        y_points_d2 = [mother_y, daughter_2_y]
        
        # plot connecting line for daughter cell 1
        plt.plot(x_points_d1, y_points_d1, color = 'k',# linestyle = 'dotted',
         linewidth = 3, markersize = 3) # line connecting start and end point black 
        plt.plot(x_points_d1, y_points_d1, color = color_i,
         linewidth=2.5, markersize=2.5) # line connecting start and end point color
        
        # plot connecting line for daughter cell 2
        plt.plot(x_points_d2, y_points_d2, color = 'k', #linestyle = 'dotted',
         linewidth = 3, markersize = 3) # line connecting start and end point black 
        plt.plot(x_points_d2, y_points_d2, color = color_i,
         linewidth=2.5, markersize=2.5) # line connecting start and end point color
        
        # start point daughter cell 1
        plt.plot(cell_list[i_sel][4][0],cell_list[i_sel][5][0], color = 'k',
          marker='o', linewidth=20, markersize=23) # end point black edge
        plt.plot(cell_list[i_sel][4][0],cell_list[i_sel][5][0], color = color_i,
          marker='o', linewidth=20, markersize=21) # end point colored
        plt.plot(cell_list[i_sel][4][0],cell_list[i_sel][5][0], color = 'k',
          marker='o', linewidth=20, markersize=16) #14 end point black
        plt.plot(cell_list[i_sel][4][0],cell_list[i_sel][5][0], color = 'w', marker=cell_marker_i,
         linewidth=10, markersize=10*multiplying_factor) # end point number in white
        
        # start point daughter cell 2
        plt.plot(cell_list[i_sel][6][0],cell_list[i_sel][7][0], color = 'k',
          marker='o', linewidth=20, markersize=23) # end point black edge
        plt.plot(cell_list[i_sel][6][0],cell_list[i_sel][7][0], color = color_i,
          marker='o', linewidth=20, markersize=21) # end point colored
        plt.plot(cell_list[i_sel][6][0],cell_list[i_sel][7][0], color = 'k',
          marker='o', linewidth=20, markersize=16) # end point black
        plt.plot(cell_list[i_sel][6][0],cell_list[i_sel][7][0], color = 'w', marker=cell_marker_i,
         linewidth=20, markersize=10*multiplying_factor) # end point number
    
    # start points non-dividing cells and mother cells
    plt.plot(cell_list[i_sel][0][0],cell_list[i_sel][1][0], color = 'k', marker='o',
     linewidth=2, markersize=23) # start point black edge
    plt.plot(cell_list[i_sel][0][0],cell_list[i_sel][1][0], color = color_i, marker='o',
     linewidth=2, markersize=21) # start point colored
    plt.plot(cell_list[i_sel][0][0],cell_list[i_sel][1][0], color = 'w', marker='o',
     linewidth=2, markersize=16) # start point white
    plt.plot(cell_list[i_sel][0][0],cell_list[i_sel][1][0], color = 'k', marker=cell_marker_i,
     linewidth=2, markersize=10*multiplying_factor) # start point number
    
    # end points
    end_points = 1
    if end_points == 1:
        if cell_list[i_sel][2] == 'Dividing cell':
            # daughter cell 1
            plt.plot(cell_list[i_sel][4][-1],cell_list[i_sel][5][-1], color = 'k',
              marker='D', linewidth=20, markersize=22) # end point black edge
            plt.plot(cell_list[i_sel][4][-1],cell_list[i_sel][5][-1], color = color_i,
              marker='D', linewidth=20, markersize=20) # end point colored
            plt.plot(cell_list[i_sel][4][-1],cell_list[i_sel][5][-1], color = 'w',
              marker='D', linewidth=20, markersize=15) # end point white
            plt.plot(cell_list[i_sel][4][-1],cell_list[i_sel][5][-1], color = 'k', marker=cell_marker_i,
             linewidth=20, markersize=8*multiplying_factor) # end point number
        
            # daughter cell 2
            plt.plot(cell_list[i_sel][6][-1],cell_list[i_sel][7][-1], color = 'k',
              marker='D', linewidth=20, markersize=22) # end point black edge
            plt.plot(cell_list[i_sel][6][-1],cell_list[i_sel][7][-1], color = color_i,
              marker='D', linewidth=20, markersize=20) # end point colored
            plt.plot(cell_list[i_sel][6][-1],cell_list[i_sel][7][-1], color = 'w',
              marker='D', linewidth=20, markersize=15) #13 end point white
            plt.plot(cell_list[i_sel][6][-1],cell_list[i_sel][7][-1], color = 'k', marker=cell_marker_i,
             linewidth=20, markersize=10*multiplying_factor) # end point number
        
        else: # no dividing cell
            plt.plot(cell_list[i_sel][0][-1],cell_list[i_sel][1][-1], color = 'k',
              marker='D', linewidth=20, markersize=22) # end point black edge
            plt.plot(cell_list[i_sel][0][-1],cell_list[i_sel][1][-1], color = color_i,
              marker='D', linewidth=20, markersize=20) # end point colored
            plt.plot(cell_list[i_sel][0][-1],cell_list[i_sel][1][-1], color = 'w',
              marker='D', linewidth=20, markersize=15) # end point white
            plt.plot(cell_list[i_sel][0][-1],cell_list[i_sel][1][-1], color = 'k', marker=cell_marker_i,
             linewidth=20, markersize=10*multiplying_factor) # end point number

#get current axes
ax = plt.gca()

#hide x-axis
ax.get_xaxis().set_visible(False)

#hide y-axis 
ax.get_yaxis().set_visible(False)

# Save figure
#plt.savefig("Mosaic_TL_mCherry/Figure_in_vitro_tracking_mosaic_icons-t31_continuous-line_black.png", dpi=600, bbox_inches='tight',pad_inches = 0)

# Plot figure
fig = plt.show()


# CORRECTION FOR DRIFT only
# For T = 31, plot 14 cell trajectories that are corrected for the 'drift' of the center's gastruloid (center of ellipse) over time
# With drift we mean: how the center of the gastruloid - estimated by the center of the fitted ellipse - changes in position over time w.r.t. to the center of the gastruloid at the final time step.
%matplotlib inline
plt.rcParams['figure.figsize'] = [30, 20]

imgplot= mpimg.imread('Mosaic_TL_mCherry/MAX_prediction-t31-200epochsRGB-SB100um.png')
plt.imshow(imgplot)#, extent=[0, 300, 0, 300])

cmap = [#'#000000', #black
        '#e69d00', #orange
        '#f0e442', #yellow
        '#009e73', #bluish green
        '#56b4e9', #sky blue
        '#cc79a7', #reddish purple
        '#d55e00', #vermillion (reddish)
        '#0072b2', #blue
        #'#ffffff', #white
        '#7f7f7f'] #grey
cell_marker = ['$1$','$2$','$3$','$4$','$5$','$6$','$7$','$8$','$9$','$10$','$11$','$12$','$13$','$14$']
n = 14

# Calculate drift
x_corr_drift_iv = [None] * t_steps
y_corr_drift_iv = [None] * t_steps
for j in range(0,t_steps):
    x_corr_drift_iv[j] = xc_iv[-1]-xc_iv[j]
    y_corr_drift_iv[j] = yc_iv[-1]-yc_iv[j]

# Plot trajectories
for i in range(n):
    if i >= 9:
        multiplying_factor = 1.5 # use a larger icon for cells that start to count from 10
    else:
        multiplying_factor = 1
       
    i_sel = i
    cell_marker_i = cell_marker[i_sel]
    if i < 8:
        color_i = cmap[i]
    if i >= 8 and i < 16:
        color_i = cmap[i-8]

    #Trajectories
    # Dividing cells
    # division points; cell 3: 17, cell 5: 16, cell 7: 14, cell 9: 12
    if cell_list[i_sel][2] == 'Dividing cell':
        cell_division_point = cell_list[i_sel][3]
        
        plt.plot(cell_list[i_sel][4]+x_corr_drift_iv[-len(cell_list[i_sel][4]):], cell_list[i_sel][5]+y_corr_drift_iv[-len(cell_list[i_sel][5]):], color = 'k',
         linewidth = 3, markersize = 3) # black trajectory "daughter cell 1"
        plt.plot(cell_list[i_sel][4]+x_corr_drift_iv[-len(cell_list[i_sel][4]):], cell_list[i_sel][5]+y_corr_drift_iv[-len(cell_list[i_sel][5]):], color = color_i,
          linewidth=2.5, markersize=2.5) # colored trajectory "daughter cell 1"
        
        plt.plot(cell_list[i_sel][6]+x_corr_drift_iv[-len(cell_list[i_sel][6]):], cell_list[i_sel][7]+y_corr_drift_iv[-len(cell_list[i_sel][7]):], color = 'k',
         linewidth = 3, markersize = 3) # black trajectory "daughter cell 2"
        plt.plot(cell_list[i_sel][6]+x_corr_drift_iv[-len(cell_list[i_sel][6]):], cell_list[i_sel][7]+y_corr_drift_iv[-len(cell_list[i_sel][7]):], color = color_i,
          linewidth=2.5, markersize=2.5) # colored trajectory "daughter cell 2"
       
    # Mother cell
    # Trajectory from start to end
    plt.plot(cell_list[i_sel][0]+x_corr_drift_iv[0:len(cell_list[i_sel][0])], cell_list[i_sel][1]+y_corr_drift_iv[0:len(cell_list[i_sel][1])], color = 'k',
     linewidth = 3, markersize = 3) # line connecting start and end point black 
    plt.plot(cell_list[i_sel][0]+x_corr_drift_iv[0:len(cell_list[i_sel][0])], cell_list[i_sel][1]+y_corr_drift_iv[0:len(cell_list[i_sel][1])], color = color_i,
     linewidth=2.5, markersize=2.5) # line connecting start and end point
    
    # Dividing cells - connecting line and starting points
    if cell_list[i_sel][2] == 'Dividing cell':
        
        # Connecting line mother cell and daughter cell
        mother_x = cell_list[i_sel][0][-1]+x_corr_drift_iv[len(cell_list[i_sel][0])-1]
        mother_y = cell_list[i_sel][1][-1]+y_corr_drift_iv[len(cell_list[i_sel][1])-1]
        daughter_1_x = cell_list[i_sel][4][0]+x_corr_drift_iv[len(cell_list[i_sel][0])]
        daughter_1_y = cell_list[i_sel][5][0]+y_corr_drift_iv[len(cell_list[i_sel][1])]
        daughter_2_x = cell_list[i_sel][6][0]+x_corr_drift_iv[len(cell_list[i_sel][0])]
        daughter_2_y = cell_list[i_sel][7][0]+y_corr_drift_iv[len(cell_list[i_sel][1])]
        x_points_d1 = [mother_x, daughter_1_x]
        x_points_d2 = [mother_x, daughter_2_x]
        y_points_d1 = [mother_y, daughter_1_y]
        y_points_d2 = [mother_y, daughter_2_y]
        
        plt.plot(x_points_d1, y_points_d1, color = 'k',# linestyle = 'dotted',
         linewidth = 3, markersize = 3) # line connecting start and end point black 
        plt.plot(x_points_d1, y_points_d1, color = color_i,
         linewidth=2.5, markersize=2.5) # line connecting start and end point
        
        plt.plot(x_points_d2, y_points_d2, color = 'k', #linestyle = 'dotted',
         linewidth = 3, markersize = 3) # line connecting start and end point black 
        plt.plot(x_points_d2, y_points_d2, color = color_i,
         linewidth=2.5, markersize=2.5) # line connecting start and end point
        
        # Start point daughter cell 1
        plt.plot(cell_list[i_sel][4][0]+x_corr_drift_iv[len(cell_list[i_sel][0])],cell_list[i_sel][5][0]+y_corr_drift_iv[len(cell_list[i_sel][1])], color = 'k',
          marker='o', linewidth=20, markersize=23) # end point black edge
        plt.plot(cell_list[i_sel][4][0]+x_corr_drift_iv[len(cell_list[i_sel][0])],cell_list[i_sel][5][0]+y_corr_drift_iv[len(cell_list[i_sel][1])], color = color_i,
          marker='o', linewidth=20, markersize=21) # end point colored
        plt.plot(cell_list[i_sel][4][0]+x_corr_drift_iv[len(cell_list[i_sel][0])],cell_list[i_sel][5][0]+y_corr_drift_iv[len(cell_list[i_sel][1])], color = 'k',
          marker='o', linewidth=20, markersize=16) #14 end point black
        plt.plot(cell_list[i_sel][4][0]+x_corr_drift_iv[len(cell_list[i_sel][0])],cell_list[i_sel][5][0]+y_corr_drift_iv[len(cell_list[i_sel][1])], color = 'w', marker=cell_marker_i,
         linewidth=10, markersize=10*multiplying_factor) # end point number in white
        
        # Start point daughter cell 2
        plt.plot(cell_list[i_sel][6][0]+x_corr_drift_iv[len(cell_list[i_sel][0])],cell_list[i_sel][7][0]+y_corr_drift_iv[len(cell_list[i_sel][1])], color = 'k',
          marker='o', linewidth=20, markersize=23) # end point black edge
        plt.plot(cell_list[i_sel][6][0]+x_corr_drift_iv[len(cell_list[i_sel][0])],cell_list[i_sel][7][0]+y_corr_drift_iv[len(cell_list[i_sel][1])], color = color_i,
          marker='o', linewidth=20, markersize=21) # end point colored
        plt.plot(cell_list[i_sel][6][0]+x_corr_drift_iv[len(cell_list[i_sel][0])],cell_list[i_sel][7][0]+y_corr_drift_iv[len(cell_list[i_sel][1])], color = 'k',
          marker='o', linewidth=20, markersize=16) # end point black
        plt.plot(cell_list[i_sel][6][0]+x_corr_drift_iv[len(cell_list[i_sel][0])],cell_list[i_sel][7][0]+y_corr_drift_iv[len(cell_list[i_sel][1])], color = 'w', marker=cell_marker_i,
         linewidth=20, markersize=10*multiplying_factor) # end point number
   
    # Starting points non-dividing cells    
    plt.plot(cell_list[i_sel][0][0]+x_corr_drift_iv[0],cell_list[i_sel][1][0]+y_corr_drift_iv[0], color = 'k', marker='o',
     linewidth=2, markersize=23) # start point black edge
    plt.plot(cell_list[i_sel][0][0]+x_corr_drift_iv[0],cell_list[i_sel][1][0]+y_corr_drift_iv[0], color = color_i, marker='o',
     linewidth=2, markersize=21) # start point colored
    plt.plot(cell_list[i_sel][0][0]+x_corr_drift_iv[0],cell_list[i_sel][1][0]+y_corr_drift_iv[0], color = 'w', marker='o',
     linewidth=2, markersize=16) # start point white
    plt.plot(cell_list[i_sel][0][0]+x_corr_drift_iv[0],cell_list[i_sel][1][0]+y_corr_drift_iv[0], color = 'k', marker=cell_marker_i,
     linewidth=2, markersize=10*multiplying_factor) # start point number
    
    # End points (all cells)
    end_points = 1
    if end_points == 1:
        if cell_list[i_sel][2] == 'Dividing cell':
            # daughter cell 1
            plt.plot(cell_list[i_sel][4][-1],cell_list[i_sel][5][-1], color = 'k',
              marker='D', linewidth=20, markersize=22) # end point black edge
            plt.plot(cell_list[i_sel][4][-1],cell_list[i_sel][5][-1], color = color_i,
              marker='D', linewidth=20, markersize=20) # end point colored
            plt.plot(cell_list[i_sel][4][-1],cell_list[i_sel][5][-1], color = 'w',
              marker='D', linewidth=20, markersize=15) # end point white
            plt.plot(cell_list[i_sel][4][-1],cell_list[i_sel][5][-1], color = 'k', marker=cell_marker_i,
             linewidth=20, markersize=8*multiplying_factor) # end point number
        
            # daughter cell 2
            plt.plot(cell_list[i_sel][6][-1],cell_list[i_sel][7][-1], color = 'k',
              marker='D', linewidth=20, markersize=22) # end point black edge
            plt.plot(cell_list[i_sel][6][-1],cell_list[i_sel][7][-1], color = color_i,
              marker='D', linewidth=20, markersize=20) # end point colored
            plt.plot(cell_list[i_sel][6][-1],cell_list[i_sel][7][-1], color = 'w',
              marker='D', linewidth=20, markersize=15) #13 end point white
            plt.plot(cell_list[i_sel][6][-1],cell_list[i_sel][7][-1], color = 'k', marker=cell_marker_i,
             linewidth=20, markersize=10*multiplying_factor) # end point number
        
        else: # no dividing cell
            plt.plot(cell_list[i_sel][0][-1],cell_list[i_sel][1][-1], color = 'k',
              marker='D', linewidth=20, markersize=22) # end point black edge
            plt.plot(cell_list[i_sel][0][-1],cell_list[i_sel][1][-1], color = color_i,
              marker='D', linewidth=20, markersize=20) # end point colored
            plt.plot(cell_list[i_sel][0][-1],cell_list[i_sel][1][-1], color = 'w',
              marker='D', linewidth=20, markersize=15) # end point white
            plt.plot(cell_list[i_sel][0][-1],cell_list[i_sel][1][-1], color = 'k', marker=cell_marker_i,
             linewidth=20, markersize=10*multiplying_factor) # end point number

#get current axes
ax = plt.gca()

#hide x-axis
ax.get_xaxis().set_visible(False)

#hide y-axis 
ax.get_yaxis().set_visible(False)

# Save figure
#plt.savefig("Mosaic_TL_mCherry/Figure_in_vitro_tracking_mosaic_icons-t31_continuous-line_black.png", dpi=600, bbox_inches='tight',pad_inches = 0)

# Plot figure
fig = plt.show()


# CORRECTION FOR DRIFT AND ROTATION
# For T = 31, plot 14 cells trajectories that are corrected for the changing rotation of the gastruloid and the 'drift' of the center's gastruloid (center of ellipse) over time.
# With drift we mean: how the center of the gastruloid - estimated by the center of the fitted ellipse - changes in position over time w.r.t. to the center of the gastruloid at the final time step.
# With rotation we mean: how the rotation of the gastruloid changes w.r.t. the gastruloid's orientation at the final time step. To correct for the rotation we use the rotated cell centers and
# rotate them again to the orientation of the gastruloid at the final time step.

%matplotlib inline
plt.rcParams['figure.figsize'] = [30, 20]

imgplot= mpimg.imread('Mosaic_TL_mCherry/MAX_prediction-t31-200epochsRGB-SB100um.png')
plt.imshow(imgplot)#, extent=[0, 300, 0, 300])

# colors per cell
cmap = [#'#000000', #black
        '#e69d00', #orange # 1
        '#f0e442', #yellow # 2
        '#009e73', #bluish green # 3a
        '#009e73', #bluish green # 3b
        '#56b4e9', #sky blue # 4
        '#cc79a7', #reddish purple # 5a
        '#cc79a7', #reddish purple # 5b
        '#d55e00', #vermillion (reddish) # 6
        '#0072b2', #blue # 7a
        '#0072b2', #blue # 7b
        #'#ffffff', #white
        '#7f7f7f', #grey # 8
        '#e69d00', #orange # 9a
        '#e69d00', #orange # 9b
        '#f0e442', #yellow # 10
        '#009e73', #bluish green # 11
        '#56b4e9', #sky blue # 12
        '#cc79a7', #reddish purple # 13
        '#d55e00'] #vermillion (reddish) # 14

cell_marker = ['$1$','$2$','$3$','$3$','$4$','$5$','$5$','$6$','$7$','$7$','$8$','$9$','$9$','$10$','$11$','$12$','$13$','$14$']

# Rotate the rotated cell centers to the orientation of the gastruloid at the final time step
no_cells = len(rotated_cell_centers_iv)
rotated_cell_centers_rot_iv = np.empty((no_cells,t_steps,4))
for i in range(no_cells): # cells

    for j in range(t_steps): # time
        # O = position of cell
        O_x_rot_iv = np.cos(2*np.pi - theta_corr_iv[-1]) * (rotated_cell_centers_iv[i,j,0]-xc_iv[j]) + np.sin(2*np.pi - theta_corr_iv[-1]) * (rotated_cell_centers_iv[i,j,1]-yc_iv[j])
        rotated_cell_centers_rot_iv[i,j,0] = O_x_rot_iv + xc_iv[j] # rotated x-coordinate of cell
        rotated_cell_centers_rot_iv[i,j,2] = O_x_rot_iv # x component of vector from center of ellipse to O

        O_y_rot_iv = np.cos(2*np.pi - theta_corr_iv[-1]) * (rotated_cell_centers_iv[i,j,1]-yc_iv[j]) - np.sin(2*np.pi - theta_corr_iv[-1]) * (rotated_cell_centers_iv[i,j,0]-xc_iv[j])
        rotated_cell_centers_rot_iv[i,j,1] = O_y_rot_iv + yc_iv[j] # rotated y-coordinate of cell
        rotated_cell_centers_rot_iv[i,j,3] = O_y_rot_iv # y component of vector from center of ellipse to O

# Calculate drift
x_corr_drift_iv = [None] * t_steps
y_corr_drift_iv = [None] * t_steps

for j in range(0,t_steps):
    x_corr_drift_iv[j] = xc_iv[-1]-xc_iv[j]
    y_corr_drift_iv[j] = yc_iv[-1]-yc_iv[j]
    
# Plot trajectories corrected for DRIFT (add drift) + corrected for ROTATION (plot rotated-rotated cell centers)
for i in range(no_cells):
    
    # icon number size and color
    if i > 12:
        multiplying_factor = 1.5 # use a larger icon for cells that start to count from 10
    else:
        multiplying_factor = 1
    cell_marker_i = cell_marker[i]
    color_i = cmap[i]
    
    # Plot trajectories (rotated-rotated) and add drift
    plt.plot(rotated_cell_centers_rot_iv[i,:,0]+x_corr_drift_iv, rotated_cell_centers_rot_iv[i,:,1]+y_corr_drift_iv, color = 'k', #linestyle = 'dotted',
     linewidth = 3, markersize = 3) # line connecting start and end point black 
    plt.plot(rotated_cell_centers_rot_iv[i,:,0]+x_corr_drift_iv, rotated_cell_centers_rot_iv[i,:,1]+y_corr_drift_iv, color = color_i,
     linewidth=2.5, markersize=2.5) # line connecting start and end point   
    
# Plot starting points
for i in range(no_cells):
    
    # icon number size and color
    if i > 12:
        multiplying_factor = 1.5 # use a larger icon for cells that start to count from 10
    else:
        multiplying_factor = 1
    cell_marker_i = cell_marker[i]
    color_i = cmap[i]
    
    # For all cells, plot start icon
    plt.plot(rotated_cell_centers_rot_iv[i,0,0]+x_corr_drift_iv[0],rotated_cell_centers_rot_iv[i,0,1]+y_corr_drift_iv[0], color = 'k',
     marker='o', linewidth=20, markersize=23) # end point black edge
    plt.plot(rotated_cell_centers_rot_iv[i,0,0]+x_corr_drift_iv[0],rotated_cell_centers_rot_iv[i,0,1]+y_corr_drift_iv[0], color = color_i,
     marker='o', linewidth=20, markersize=21) # end point colored
    plt.plot(rotated_cell_centers_rot_iv[i,0,0]+x_corr_drift_iv[0],rotated_cell_centers_rot_iv[i,0,1]+y_corr_drift_iv[0], color = 'w',
     marker='o', linewidth=20, markersize=16) #14 point white
    plt.plot(rotated_cell_centers_rot_iv[i,0,0]+x_corr_drift_iv[0],rotated_cell_centers_rot_iv[i,0,1]+y_corr_drift_iv[0], color = 'k', marker=cell_marker_i,
     linewidth=10, markersize=10*multiplying_factor) # number in black
  
    # For dividing cells, only plot dividing cell (end point mother cell)
    if i == 2: # cell 3
        sel_cell = i
        div_no = cell_list[sel_cell][3] - 2# 17 > 16
    if i == 5: # cell 5
        sel_cell = i - 1
        div_no = cell_list[sel_cell][3] - 2
    if i == 8: # cell 7
        sel_cell = i - 2
        div_no = cell_list[sel_cell][3] - 2
    if i == 11: # cell 9
        sel_cell = i - 3
        div_no = cell_list[sel_cell][3] - 2

    if i == 2 or i == 5 or i == 8 or i == 11:
        #End point mother cell
        plt.plot(rotated_cell_centers_rot_iv[i,div_no,0]+x_corr_drift_iv[div_no],rotated_cell_centers_rot_iv[i,div_no,1]+y_corr_drift_iv[div_no], color = 'k',
         marker='o', linewidth=20, markersize=23) # end point black edge
        plt.plot(rotated_cell_centers_rot_iv[i,div_no,0]+x_corr_drift_iv[div_no],rotated_cell_centers_rot_iv[i,div_no,1]+y_corr_drift_iv[div_no], color = color_i,
         marker='o', linewidth=20, markersize=21) # end point colored
        plt.plot(rotated_cell_centers_rot_iv[i,div_no,0]+x_corr_drift_iv[div_no],rotated_cell_centers_rot_iv[i,div_no,1]+y_corr_drift_iv[div_no], color = 'k',
         marker='o', linewidth=20, markersize=16) #14 end point black
        plt.plot(rotated_cell_centers_rot_iv[i,div_no,0]+x_corr_drift_iv[div_no],rotated_cell_centers_rot_iv[i,div_no,1]+y_corr_drift_iv[div_no], color = 'w', marker=cell_marker_i,
         linewidth=10, markersize=10*multiplying_factor) # end point number in white
   
    # For all cells, plot end icon
    plt.plot(rotated_cell_centers_rot_iv[i,-1,0]+x_corr_drift_iv[-1],rotated_cell_centers_rot_iv[i,-1,1]+y_corr_drift_iv[-1], color = 'k',
     marker='D', linewidth=20, markersize=23) # end point black edge
    plt.plot(rotated_cell_centers_rot_iv[i,-1,0]+x_corr_drift_iv[-1],rotated_cell_centers_rot_iv[i,-1,1]+y_corr_drift_iv[-1], color = color_i,
     marker='D', linewidth=20, markersize=21) # end point colored
    plt.plot(rotated_cell_centers_rot_iv[i,-1,0]+x_corr_drift_iv[-1],rotated_cell_centers_rot_iv[i,-1,1]+y_corr_drift_iv[-1], color = 'w',
     marker='D', linewidth=20, markersize=16) #14 point white
    plt.plot(rotated_cell_centers_rot_iv[i,-1,0]+x_corr_drift_iv[-1],rotated_cell_centers_rot_iv[i,-1,1]+y_corr_drift_iv[-1], color = 'k', marker=cell_marker_i,
     linewidth=10, markersize=10*multiplying_factor) # number in black
            
#get current axes
ax = plt.gca()

#hide x-axis
ax.get_xaxis().set_visible(False)

#hide y-axis 
ax.get_yaxis().set_visible(False)

# Save figure
#plt.savefig("Mosaic_TL_mCherry/Figure_in_vitro_tracking_mosaic_t31_corrected_drift_rotation.png", dpi=600, bbox_inches='tight',pad_inches = 0)

# Plot figure
fig = plt.show()


# T = 1, 14 cells starting points + yellow connecting lines
plt.rcParams['figure.figsize'] = [30, 20]

imgplot= mpimg.imread('Mosaic_TL_mCherry/MAX_prediction-t1-200epochsRGB_SB100um.png')
plt.imshow(imgplot)#, extent=[0, 300, 0, 300])


cmap = [#'#000000', #black
        '#e69d00', #orange
        '#f0e442', #yellow
        '#009e73', #bluish green
        '#56b4e9', #sky blue
        '#cc79a7', #reddish purple
        '#d55e00', #vermillion (reddish)
        '#0072b2', #blue
        #'#ffffff', #white
        '#7f7f7f'] #grey
cell_marker = ['$1$','$2$','$3$','$4$','$5$','$6$','$7$','$8$','$9$','$10$','$11$','$12$','$13$','$14$']
n = 14

# extra lines

#border
x_points_lines = [cell_list[6][0][0],cell_list[2][0][0],cell_list[1][0][0],cell_list[0][0][0],
                  cell_list[5][0][0],cell_list[11][0][0],cell_list[13][0][0],cell_list[12][0][0],cell_list[9][0][0],cell_list[6][4][0]]
y_points_lines = [cell_list[6][1][0],cell_list[2][1][0],cell_list[1][1][0],cell_list[0][1][0],
                  cell_list[5][1][0],cell_list[11][1][0],cell_list[13][1][0],cell_list[12][1][0],cell_list[9][1][0],cell_list[6][1][0]]
#in between
x_points_lines_mid_1 = [cell_list[2][0][0],cell_list[3][0][0],cell_list[5][0][0]]
y_points_lines_mid_1 = [cell_list[2][1][0],cell_list[3][1][0],cell_list[5][1][0]]
x_points_lines_mid_2 = [cell_list[6][0][0],cell_list[7][0][0],cell_list[5][0][0]]
y_points_lines_mid_2 = [cell_list[6][1][0],cell_list[7][1][0],cell_list[5][1][0]]
x_points_lines_mid_3 = [cell_list[9][0][0],cell_list[10][0][0],cell_list[11][0][0]]
y_points_lines_mid_3 = [cell_list[9][1][0],cell_list[10][1][0],cell_list[11][1][0]]

plt.plot(x_points_lines, y_points_lines,'y',linewidth = 5)
plt.plot(x_points_lines_mid_1, y_points_lines_mid_1,'y',linewidth = 5)
plt.plot(x_points_lines_mid_2, y_points_lines_mid_2,'y',linewidth = 5)
plt.plot(x_points_lines_mid_3, y_points_lines_mid_3,'y',linewidth = 5)

for i in range(n):
    if i >= 9:
        multiplying_factor = 1.5 # use a larger icon for cells that start to count from 10
    else:
        multiplying_factor = 1
       
    i_sel = i
    cell_marker_i = cell_marker[i_sel]
    if i < 8:
        color_i = cmap[i]
    if i >= 8 and i < 16:
        color_i = cmap[i-8]

    # start points   
    plt.plot(cell_list[i_sel][0][0],cell_list[i_sel][1][0], color = 'k', marker='o',
     linewidth=2, markersize=23) # start point black edge
    plt.plot(cell_list[i_sel][0][0],cell_list[i_sel][1][0], color = color_i, marker='o',
     linewidth=2, markersize=21) # start point colored
    plt.plot(cell_list[i_sel][0][0],cell_list[i_sel][1][0], color = 'w', marker='o',
     linewidth=2, markersize=16) # start point white
    plt.plot(cell_list[i_sel][0][0],cell_list[i_sel][1][0], color = 'k', marker=cell_marker_i,
     linewidth=2, markersize=10*multiplying_factor) # start point number
     
#get current axes
ax = plt.gca()

#hide x-axis
ax.get_xaxis().set_visible(False)

#hide y-axis 
ax.get_yaxis().set_visible(False)

# Save figure
#plt.savefig("Mosaic_TL_mCherry/Figure_in_vitro_tracking_mosaic_icons_lines_t1_more_white.png", dpi=600, bbox_inches='tight',pad_inches = 0)

# Plot figure
fig = plt.show()


# T = 31, 14 cells trajectories
%matplotlib inline
plt.rcParams['figure.figsize'] = [30, 20]

imgplot= mpimg.imread('Mosaic_TL_mCherry/MAX_prediction-t31-200epochsRGB-SB100um.png')
plt.imshow(imgplot)#, extent=[0, 300, 0, 300])


cmap = [#'#000000', #black
        '#e69d00', #orange
        '#f0e442', #yellow
        '#009e73', #bluish green
        '#56b4e9', #sky blue
        '#cc79a7', #reddish purple
        '#d55e00', #vermillion (reddish)
        '#0072b2', #blue
        #'#ffffff', #white
        '#7f7f7f'] #grey
cell_marker = ['$1$','$2$','$3$','$4$','$5$','$6$','$7$','$8$','$9$','$10$','$11$','$12$','$13$','$14$']
n = 14

# extra lines

#border
x_points_lines = [cell_list[6][4][-1],cell_list[2][4][-1],cell_list[1][0][-1],cell_list[0][0][-1],
                  cell_list[5][0][-1],cell_list[11][0][-1],cell_list[13][0][-1],cell_list[12][0][-1],cell_list[9][0][-1],cell_list[6][4][-1]]
y_points_lines = [cell_list[6][5][-1],cell_list[2][5][-1],cell_list[1][1][-1],cell_list[0][1][-1],
                  cell_list[5][1][-1],cell_list[11][1][-1],cell_list[13][1][-1],cell_list[12][1][-1],cell_list[9][1][-1],cell_list[6][5][-1]]
#in between
x_points_lines_mid_1 = [cell_list[2][4][-1],cell_list[3][0][-1],cell_list[5][0][-1]]
y_points_lines_mid_1 = [cell_list[2][5][-1],cell_list[3][1][-1],cell_list[5][1][-1]]
x_points_lines_mid_2 = [cell_list[6][4][-1],cell_list[7][0][-1],cell_list[5][0][-1]]
y_points_lines_mid_2 = [cell_list[6][5][-1],cell_list[7][1][-1],cell_list[5][1][-1]]
x_points_lines_mid_3 = [cell_list[9][0][-1],cell_list[10][0][-1],cell_list[11][0][-1]]
y_points_lines_mid_3 = [cell_list[9][1][-1],cell_list[10][1][-1],cell_list[11][1][-1]]

plt.plot(x_points_lines, y_points_lines,'y',linewidth = 5)
plt.plot(x_points_lines_mid_1, y_points_lines_mid_1,'y',linewidth = 5)
plt.plot(x_points_lines_mid_2, y_points_lines_mid_2,'y',linewidth = 5)
plt.plot(x_points_lines_mid_3, y_points_lines_mid_3,'y',linewidth = 5)

for i in range(n):
    if i >= 9:
        multiplying_factor = 1.5 # use a larger icon for cells that start to count from 10
    else:
        multiplying_factor = 1
       
    i_sel = i
    cell_marker_i = cell_marker[i_sel]
    if i < 8:
        color_i = cmap[i]
    if i >= 8 and i < 16:
        color_i = cmap[i-8]

    # end points
    end_points = 1
    if end_points == 1:
        if cell_list[i_sel][2] == 'Dividing cell':
            # daughter cell 1
            plt.plot(cell_list[i_sel][4][-1],cell_list[i_sel][5][-1], color = 'k',
              marker='D', linewidth=20, markersize=22) # end point black edge
            plt.plot(cell_list[i_sel][4][-1],cell_list[i_sel][5][-1], color = color_i,
              marker='D', linewidth=20, markersize=20) # end point colored
            plt.plot(cell_list[i_sel][4][-1],cell_list[i_sel][5][-1], color = 'w',
              marker='D', linewidth=20, markersize=15) # end point white
            plt.plot(cell_list[i_sel][4][-1],cell_list[i_sel][5][-1], color = 'k', marker=cell_marker_i,
             linewidth=20, markersize=8*multiplying_factor) # end point number
        
            # daughter cell 2
            plt.plot(cell_list[i_sel][6][-1],cell_list[i_sel][7][-1], color = 'k',
              marker='D', linewidth=20, markersize=22) # end point black edge
            plt.plot(cell_list[i_sel][6][-1],cell_list[i_sel][7][-1], color = color_i,
              marker='D', linewidth=20, markersize=20) # end point colored
            plt.plot(cell_list[i_sel][6][-1],cell_list[i_sel][7][-1], color = 'w',
              marker='D', linewidth=20, markersize=15) # end point white
            plt.plot(cell_list[i_sel][6][-1],cell_list[i_sel][7][-1], color = 'k', marker=cell_marker_i,
             linewidth=20, markersize=10*multiplying_factor) # end point number
        
        else: # no dividing cell
            plt.plot(cell_list[i_sel][0][-1],cell_list[i_sel][1][-1], color = 'k',
              marker='D', linewidth=20, markersize=22) # end point black edge
            plt.plot(cell_list[i_sel][0][-1],cell_list[i_sel][1][-1], color = color_i,
              marker='D', linewidth=20, markersize=20) # end point colored
            plt.plot(cell_list[i_sel][0][-1],cell_list[i_sel][1][-1], color = 'w',
              marker='D', linewidth=20, markersize=15) # end point white
            plt.plot(cell_list[i_sel][0][-1],cell_list[i_sel][1][-1], color = 'k', marker=cell_marker_i,
             linewidth=20, markersize=10*multiplying_factor) # end point number
            
#get current axes
ax = plt.gca()

#hide x-axis
ax.get_xaxis().set_visible(False)

#hide y-axis 
ax.get_yaxis().set_visible(False)

# Save figure
#plt.savefig("Figure_in_vitro_tracking_mosaic_icons_lines_t31_endpointsonly.png", dpi=600, bbox_inches='tight',pad_inches = 0)

# Plot figure
fig = plt.show()


t_end_iv = 30 # t = 31
t_start_iv = 0 # t = 1

# Determine the length of the gastruloid at t_end_iv = 2 * long arm of ellipse fit. The long arm should be stored in the
# parameter 'b', and the short arm in parameter 'a', but the fitting tool sometimes stored them incorrectly.
# Therefore, check the value of a and b first to determine which parameter corresponds to the long arm.
# In our in vitro example, the long arm of the fitted ellipses should be on the y-axis. If a and b are incorrectly switched,
# the values stored in the x coordinates correspond to the long arm.

if a_iv[t_end_iv] > b_iv[t_end_iv]:
    
    length_gastruloid_iv = 2*a_iv[t_end_iv] # used to normalize the distance the cells moved along the long axis
    print('a > b, so we use the rotated x-coordinates to calculate the distance cells moved along the long axis of the gastruloid')
     
    end_pos_on_long_axis_iv = rotated_cell_centers_iv[:,t_end_iv,0] # rotated x-coordinate at t = t_end_iv
    end_pos_on_long_axis_iv_norm = (end_pos_on_long_axis_iv - xc_iv[t_end_iv])/(length_gastruloid_iv) # (rotated x-coordinate at t_end_iv - x_center_of_ellipse at t_end_iv) divided by length gastruloid
    moved_end_iv_norm = ((end_pos_on_long_axis_iv - rotated_cell_centers_iv[:,t_start_iv,0]) - (xc_iv[t_end_iv] - xc_iv[t_start_iv]) ) / length_gastruloid_iv # to normalize, divide by length of gastruloid

else: 
    length_gastruloid_iv = 2*b_iv[t_end_iv] # used to normalize the distance the cells moved along the long axis
    print('a < b, so we use the rotated y-coordinates to calculate the distance cells moved along the long axis of the gastruloid')
      
    end_pos_on_long_axis_iv = rotated_cell_centers_iv[:,t_end_iv,1] # rotated y-coordinate at t = t_end_iv
    end_pos_on_long_axis_iv_norm = (end_pos_on_long_axis_iv - yc_iv[t_end_iv])/(length_gastruloid_iv) # (rotated y-coordinate at t_end_iv - y_center_of_ellipse at t_end_iv) divided by length gastruloid
    moved_end_iv_norm = ((end_pos_on_long_axis_iv - rotated_cell_centers_iv[:,t_start_iv,1]) - (yc_iv[t_end_iv] - yc_iv[t_start_iv]) ) / length_gastruloid_iv # to normalize, divide by length of gastruloid
    
# -------- FIGURE ----------
# Figure parameters    
%matplotlib inline
font = {'size'   : 12}
matplotlib.rc('font', **font)
matplotlib.rcParams['font.sans-serif'] = "Arial"
matplotlib.rcParams['font.family'] = "sans-serif" 

plt.gcf().clear()
fig = plt.figure(1)
ax = fig.add_subplot(111)

plt.rcParams['figure.figsize'] = [10, 10]
plt.rcParams['xtick.labelsize'] = 'x-large'
plt.rcParams['ytick.labelsize'] = 'x-large'  
plt.xlabel('End position projected on long axis', fontsize = 16)
plt.ylabel('Movement along long axis', fontsize = 16)
plt.xlim(-0.55, 0.55)
plt.ylim(-0.51, 0.51)
plt.xticks([-0.5, -0.25, 0, 0.25, 0.5],['-L/2','-L/4','0','L/4','L/2'])
plt.yticks([-0.5, -0.25, 0, 0.25, 0.5],['-L/2','-L/4','0','L/4','L/2'])

# In vitro time lapse data
ax.scatter(end_pos_on_long_axis_iv_norm, moved_end_iv_norm, color = '#cc79a7', marker = 'o', label = 'in vitro',s = 50)

handles, labels = ax.get_legend_handles_labels()
lgd = ax.legend(handles, labels, loc='upper left', bbox_to_anchor=(1,1))
fig.set_size_inches((110/5)/2.54, (76/5)/2.54)

plt.show()

a < b, so we use the rotated y-coordinates to calculate the distance cells moved along the long axis of the gastruloid


t_end_iv = 30 # t = 31
t_start_iv = 0 # t = 1

# Determine the length of the gastruloid at t_end_iv = 2 * long arm of ellipse fit. The long arm should be stored in the
# parameter 'b', and the short arm in parameter 'a', but the fitting tool sometimes stored them incorrectly.
# Therefore, check the value of a and b first to determine which parameter corresponds to the long or short arm.
# In our in vitro example, the long arm of the fitted ellipses should be on the y-axis. The short arm on the x-axis. 
# If a and b are incorrectly switched, the values stored in the x coordinates correspond to the long arm, and y coordinates correspond to the short arm.

if a_iv[t_end_iv] > b_iv[t_end_iv]:
    
    length_gastruloid_iv_short = 2*b_iv[t_end_iv] # used to normalize the distance the cells moved along the short axis
    length_gastruloid_iv_long = 2*a_iv[t_end_iv] # used to normalize the position of the cells along the long axis
    print('a > b, so we use the rotated y-coordinates to calculate the distance cells moved along the short axis of the gastruloid')
      
    end_pos_on_short_axis_iv = rotated_cell_centers_iv[:,t_end_iv,1] # rotated y-coordinate at t = t_end_iv
    end_pos_on_long_axis_iv = rotated_cell_centers_iv[:,t_end_iv,0] # rotated x-coordinate at t = t_end_iv
    end_pos_on_short_axis_iv_norm = (end_pos_on_short_axis_iv - yc_iv[t_end_iv])/(length_gastruloid_iv_short) # (rotated y-coordinate at t_end_iv - y_center_of_ellipse at t_end_iv) divided by length gastruloid
    end_pos_on_long_axis_iv_norm = (end_pos_on_long_axis_iv - xc_iv[t_end_iv])/(length_gastruloid_iv_long) # (rotated y-coordinate at t_end_iv - y_center_of_ellipse at t_end_iv) divided by length gastruloid

    moved_end_iv_short_norm = ((end_pos_on_short_axis_iv - rotated_cell_centers_iv[:,t_start_iv,1]) - (yc_iv[t_end_iv] - yc_iv[t_start_iv]) ) / length_gastruloid_iv_short # to normalize, divide by length of gastruloid
    moved_end_iv_long_norm = ((end_pos_on_long_axis_iv - rotated_cell_centers_iv[:,t_start_iv,0]) - (xc_iv[t_end_iv] - xc_iv[t_start_iv]) ) / length_gastruloid_iv_long # to normalize, divide by length of gastruloid

else: 
    length_gastruloid_iv_short = 2*a_iv[t_end_iv] # used to normalize the distance the cells moved along the short axis
    length_gastruloid_iv_long = 2*b_iv[t_end_iv] # # used to normalize the position of the cells along the long axis
    print('a < b, so we use the rotated x-coordinates to calculate the distance cells moved along the short axis of the gastruloid')
     
    end_pos_on_short_axis_iv = rotated_cell_centers_iv[:,t_end_iv,0] # rotated x-coordinate at t = t_end_iv
    end_pos_on_long_axis_iv = rotated_cell_centers_iv[:,t_end_iv,1] # rotated y-coordinate at t = t_end_iv
    end_pos_on_short_axis_iv_norm = (end_pos_on_short_axis_iv - xc_iv[t_end_iv])/(length_gastruloid_iv_short) # (rotated x-coordinate at t_end_iv - x_center_of_ellipse at t_end_iv) divided by length gastruloid
    end_pos_on_long_axis_iv_norm = (end_pos_on_long_axis_iv - yc_iv[t_end_iv])/(length_gastruloid_iv_long) # (rotated y-coordinate at t_end_iv - y_center_of_ellipse at t_end_iv) divided by length gastruloid

    moved_end_iv_short_norm = ((end_pos_on_short_axis_iv - rotated_cell_centers_iv[:,t_start_iv,0]) - (xc_iv[t_end_iv] - xc_iv[t_start_iv]) ) / length_gastruloid_iv_short # to normalize, divide by length of gastruloid
    moved_end_iv_long_norm = ((end_pos_on_long_axis_iv - rotated_cell_centers_iv[:,t_start_iv,1]) - (yc_iv[t_end_iv] - yc_iv[t_start_iv]) ) / length_gastruloid_iv_long # to normalize, divide by length of gastruloid

# -------- FIGURE ----------
# Figure parameters    
%matplotlib inline
font = {'size'   : 12}
matplotlib.rc('font', **font)
matplotlib.rcParams['font.sans-serif'] = "Arial"
matplotlib.rcParams['font.family'] = "sans-serif" 

plt.gcf().clear()
fig = plt.figure(1)
ax = fig.add_subplot(111)

plt.rcParams['figure.figsize'] = [10, 10]
plt.rcParams['xtick.labelsize'] = 'x-large'
plt.rcParams['ytick.labelsize'] = 'x-large'  
plt.xlabel('End position projected on short axis', fontsize = 16)
plt.ylabel('Movement along short axis', fontsize = 16)
plt.xlim(-1, 1)
plt.ylim(-1, 1)
plt.xticks([-1, -0.5, 0, 0.5, 1],['S','-S/2','0','S/2', 'S'])
plt.yticks([-1,-0.5, 0, 0.5, 1],['-S','-S/2','0','S/2','S'])

ax.scatter(end_pos_on_short_axis_iv_norm, moved_end_iv_short_norm, color = '#cc79a7', marker = 'o', s = 50)#label = 'in vitro',s = 50)

fig.set_size_inches((110/5)/2.54, (76/5)/2.54)

plt.show()

a < b, so we use the rotated x-coordinates to calculate the distance cells moved along the short axis of the gastruloid


# Read image
%matplotlib inline
img = cv2.imread('Mosaic_3/image0099900.png',0)
plt.imshow(img,cmap='gray')
plt.show()


# Erode the image multiple times, and subsequently dilate the image multiple times

%matplotlib inline
# Taking a matrix of size 5 as the kernel
kernel = np.ones((5, 5), np.uint8)

# The first parameter is the original image,
# kernel is the matrix with which image is
# convolved and third parameter is the number
# of iterations, which will determine how much
# you want to erode/dilate a given image.

img_erosion = cv2.erode(img, kernel, iterations=4)
img_dilation = cv2.dilate(img_erosion, kernel, iterations=4)

imagem = cv2.bitwise_not(img_erosion)

plt.imshow(imagem,cmap='gray')
plt.show()
plt.imshow(img_dilation,cmap='gray')
plt.show()
# save image
#cv2.imwrite('Mosaic_3/image0099900-binary.png', imagem)


# read the images and store them in a list in the correct order
image_path = Path('Mosaic_corrected')
#images = list(image_path.glob('*.png'))
image_list = []
for i in range(0,100000,100):
    n = i
    STRT = "Mosaic_corrected/image"
    EXTN = ".png"
    filename = STRT + str(n).zfill(7) + EXTN
    image_list.append(imageio.imread(filename))


# from the list of images, write a GIF (standard 10 fps)
imageio.mimwrite('Mosaic_corrected/animated_mosaic_corrected.gif',image_list)


def extract_outline():

    for x in range(0,1000):
        ind = list()
        cont_all = list()
        image_loc = 'Mosaic_corrected/image'+"{:05d}".format(x)+'00.png'
        file_loc = 'Mosaic_corrected/outlines/outline_'+str(x) + '.txt'
        #file1 = open(file_loc,"a")
        mask = cv2.imread(image_loc)
        mask = 255-mask
        mask = mask[1:-1,1:-1]
        mask = mask[1:-1,1
                    :-1]
        mask = mask[1:-1,1:-1]

        if mask is not None:

            imgray = cv2.cvtColor(mask,cv2.COLOR_BGR2GRAY)
            ret,thresh = cv2.threshold(imgray,200,255,0)
            contours, im2 = cv2.findContours(thresh,cv2.RETR_TREE,cv2.CHAIN_APPROX_NONE)         
            maxsizeloc = 0
            maxsize = 0
            for i in range (len(contours)):
                if contours[i].shape[0] > maxsize:
                    maxsizeloc = i
                    maxsize = contours[i].shape[0]
            cont = np.zeros((contours[maxsizeloc].shape[0], 2))
            cont[:,0] = contours[maxsizeloc][:,0,0]
            cont[:,1] = contours[maxsizeloc][:,0,1]

            plt.plot(cont[:,0], cont[:,1])
            cont_all.append(cont)
            
            #for i in range (maxsize):
                #file1.write(str(cont_all[0][i][0]) + ',' + str(cont_all[0][i][1]) + '\n')
            #file1.close()

    return(cont_all)

extract_outline()

[array([[420., 296.],
        [420., 297.],
        [419., 298.],
        ...,
        [423., 296.],
        [422., 296.],
        [421., 296.]])]


# Ellipse fit parameters
xc = np.empty(1000) # center coordinate x of ellipse
yc = np.empty(1000) # center coordinate y of ellipse
a = np.empty(1000) # length short arm of ellipse
b = np.empty(1000) # length long arm of ellipse
theta= np.empty(1000) # rotation angle of ellipse

# Figure parameters
%matplotlib inline
fig, axs = plt.subplots(figsize = (10,10))#, sharex=True, sharey=True)#, figsize=(5,15))
cmap = [#'#000000', #black
        '#e69d00', #orange
        '#56b4e9', #sky blue
        '#009e73', #bluish green
        '#f0e442', #yellow
        '#0072b2', #blue
        '#d55e00', #vermillion (reddish)
        '#cc79a7', #reddish purple
        #'#ffffff', #white
        '#7f7f7f'] #grey

for i in range(1000):
    # Load outline coordinates for MCS i
    with open("Mosaic_corrected/outline_" + str(i) + ".txt") as file_name:
        points = np.loadtxt(file_name, delimiter=",")
    a_points = np.array(points)
    x = a_points[:, 0]
    y = a_points[:, 1]
    
    #Calculate ellipse fit from outline and store parameters
    ell = EllipseModel()
    ell.estimate(a_points)
    xc[i], yc[i], a[i], b[i], theta[i] = ell.params
    
    # Write parameters to txt file
    #file_loc = 'Mosaic_corrected/ellipse_params/MCS_'+str(i)+'-x_y_theta_a_b.txt'
    #file1 = open(file_loc,"a")
    #file1.write(str(xc[i]) + ',' + str(yc[i]) + ',' + str(theta[i]) + ',' + str(a[i]) + ',' + str(b[i]) + '\n')
    
    # Plot ellipse fit for one MCS
    t_end = 999
    if i == t_end:
        axs.scatter(x, y, color = cmap[4])
        axs.set_ylim(0,1000)
        axs.set_xlim(0,1000)
        axs.scatter(xc[i], yc[i], color=cmap[0], s=50)
        ell_patch = Ellipse((xc[i], yc[i]), 2*a[i], 2*b[i], theta[i]*180/np.pi, edgecolor=cmap[0], facecolor='none', linewidth = 3)
        axs.add_patch(ell_patch)

        if a[i] < b[i]:
            print('a',a[i],' < b',b[i], 'theta = ', theta[i]*180/np.pi)

# Save figure
#plt.savefig("Mosaic_corrected/ellipse_fit_t999.png", dpi=600, pad_inches = 0)    
    
plt.show()
print('Theta at MCS 100000 = ', theta[999]*180/np.pi)

Theta at MCS 100000 =  53.37818613761889


# Plot length of long arm and short arm over MCS - use to determine range of MCS to plot
%matplotlib inline
x1000 = list(range(0,1000))
for i in x1000:
    if a[i] > b[i]: # a is long arm, b is short arm
        plt.scatter(x1000[i],a[i], s = 0.5, color = 'k')
        plt.scatter(x1000[i],b[i], s = 0.5, color = 'b')
    elif a[i] < b[i]: # b is long arm, a is short arm
        plt.scatter(x1000[i],b[i], s = 0.5, color = 'k')
        plt.scatter(x1000[i],a[i], s = 0.5, color = 'b') 
t_start = 199
line_1 = t_start
plt.plot([line_1, line_1],[0, 220])
plt.show()


# Define number of cells and the selected time points you want to plot the trajectory for
No_cells = 20 # number of cells that are being tracked
selected_points = [t_start,249,299,349,399,449,499,549,599,649,699,749,799,849,899,949,999] # 16 steps
    
# Define list of arrays for x- and y-coordinates for each tracked cell
path_x_sim = [None] * (No_cells)
path_y_sim = [None] * (No_cells)
path_x_sim_final = [None] * No_cells
path_y_sim_final = [None] * No_cells


# Store x- and y-coordinates for each tracked cell in lists 'path_x_sim_final' and 'path_y_sim_final'
for i in range(10,int(No_cells*10+10),10):
    filepath = ""
    file_dir = ['Mosaic_corrected/cell_','.txt']
    filepath = file_dir[0]+str(i)+file_dir[1]
    with open(filepath) as f:
            for line in f:
                currentline = line.split(",")
                x_line = float(currentline[0])
                y_line = float(currentline[1])
                ii = int(i/10-1)
                path_x_sim[ii] = np.append(path_x_sim[ii],x_line)
                path_y_sim[ii] = np.append(path_y_sim[ii],y_line)
    path_x_sim_final[ii] = path_x_sim[ii][1:(len(path_x_sim[ii])+1)]
    path_y_sim_final[ii] = path_y_sim[ii][1:(len(path_x_sim[ii])+1)]


# Store x- and y-coordinates of each tracked cell for only the selected time points in lists 'x_sel_final' and 'y_sel_final'
x_sel = [None] * No_cells
y_sel = [None] * No_cells
x_sel_final = [None] * No_cells
y_sel_final = [None] * No_cells

for j in range(No_cells):
    for k in range(len(selected_points)):
        x_sel[j] = np.append(x_sel[j],2*path_x_sim_final[j][selected_points[k]])
        y_sel[j] = np.append(y_sel[j],2*path_y_sim_final[j][selected_points[k]])
    x_sel_final[j] = x_sel[j][1:(len(selected_points)+1)]
    y_sel_final[j] = y_sel[j][1:(len(selected_points)+1)]


%matplotlib inline

# Select image (time point) and adjust size (optional)
imgplot_100= plt.imread('Mosaic_corrected/image0019900.png')
plt.rcParams['figure.figsize'] = [200, 100]

# Make image gray (so no dark blue cells)
R_100, G_100, B_100 = imgplot_100[:,:,0], imgplot_100[:,:,1], imgplot_100[:,:,2]
imgGray_100 = 0.1 * R_100 + 0.1 * G_100 + 0.7*B_100

cmap = [#'#000000', #black
        '#e69d00', #orange
        '#56b4e9', #sky blue
        '#009e73', #bluish green
        '#f0e442', #yellow
        '#0072b2', #blue
        '#d55e00', #vermillion (reddish)
        '#7f7f7f', #grey
        '#cc79a7']#, #reddish purple
cell_marker = ['$1$','$2$','$3$','$4$','$5$','$6$','$7$','$8$','$9$','$10$','$11$','$12$','$13$','$14$','$15$','$16$','$17$','$18$','$19$','$20$']
markers_numbered_selected = [11,10,14,17,16,18,3,5,6,7,8,9] # 12 cells

for i in range(len(markers_numbered_selected)):
    
    if i <= 5:
        multiplying_factor = 1.5 # use a larger icon for cells that start to count from 10
    else:
        multiplying_factor = 1
    
    i_sel_rot = (markers_numbered_selected[i])*10-1
    i_sel = markers_numbered_selected[i]-1
    cell_marker_i = cell_marker[i_sel]
    if i == 0:
        color_i = cmap[2]
    if i == 2:
        color_i = cmap[7]
    if i > 0 and i < 8 and i != 2:
        color_i = cmap[i]
    if i == 8:
        color_i = cmap[0]
    if i == 9:
        color_i = cmap[2]
    if i == 10:
        color_i = cmap[1]
    if i == 11:
        color_i = cmap[3] 
    
    t_start_sel = 0
    plt.plot(x_sel_final[i_sel][t_start_sel],y_sel_final[i_sel][t_start_sel], color = 'k', marker='o',
     linewidth=2, markersize=130) # start point black edge
    plt.plot(x_sel_final[i_sel][t_start_sel],y_sel_final[i_sel][t_start_sel], color = color_i, marker='o',
     linewidth=2, markersize=120) # start point colored
    plt.plot(x_sel_final[i_sel][t_start_sel],y_sel_final[i_sel][t_start_sel], color = 'w', marker='o',
     linewidth=2, markersize=80) # start point white
    plt.plot(x_sel_final[i_sel][t_start_sel],y_sel_final[i_sel][t_start_sel], color = 'k', marker=cell_marker_i,
     linewidth=20, markersize=40*multiplying_factor) # start point number

# Save figure
plt.imshow(imgGray_100, cmap='gray')
# Save figure
#plt.savefig("Mosaic_corrected/Figure_in_silico_mosaic_corrected_12_icons_t199.png", bbox_inches='tight',pad_inches = 0)
plt.show()


%matplotlib inline

# Select image (time point) and adjust size (optional)
imgplot= plt.imread('Mosaic_corrected/image0099900.png') # final time point
plt.rcParams['figure.figsize'] = [200, 100]

# Make image gray (so no dark blue cells)
R, G, B = imgplot[:,:,0], imgplot[:,:,1], imgplot[:,:,2]
imgGray = 0.1 * R + 0.1 * G + 0.7*B #0.2989 * R + 0.5870 * G + 0.1140 * B

cmap = [#'#000000', #black
        '#e69d00', #orange
        '#56b4e9', #sky blue
        '#009e73', #bluish green
        '#f0e442', #yellow
        '#0072b2', #blue
        '#d55e00', #vermillion (reddish)
        '#7f7f7f', #grey
        '#cc79a7']#, #reddish purple
cell_marker = ['$1$','$2$','$3$','$4$','$5$','$6$','$7$','$8$','$9$','$10$','$11$','$12$','$13$','$14$','$15$','$16$','$17$','$18$','$19$','$20$']
markers_numbered_selected = [11,10,14,17,16,18,3,5,6,7,8,9]

for i in range(len(markers_numbered_selected)):
    
    if i <= 5:
        multiplying_factor = 1.5 # use a larger icon for cells that start to count from 10
    else:
        multiplying_factor = 1
    
    i_sel_rot = (markers_numbered_selected[i])*10-1
    i_sel = markers_numbered_selected[i]-1
    cell_marker_i = cell_marker[i_sel]
    if i == 0:
        color_i = cmap[2]
    if i == 2:
        color_i = cmap[7]
    if i > 0 and i < 8 and i != 2:
        color_i = cmap[i]
    if i == 8:
        color_i = cmap[0]
    if i == 9:
        color_i = cmap[2]
    if i == 10:
        color_i = cmap[1]
    if i == 11:
        color_i = cmap[3] 
   
    plt.plot(x_sel_final[i_sel][-1],y_sel_final[i_sel][-1], color = 'k',
      marker='D', linewidth=20, markersize=130) # end point black edge
    plt.plot(x_sel_final[i_sel][-1],y_sel_final[i_sel][-1], color = color_i,
      marker='D', linewidth=20, markersize=120) # end point colored
    plt.plot(x_sel_final[i_sel][-1],y_sel_final[i_sel][-1], color = 'w',
      marker='D', linewidth=20, markersize=70) # end point white
    plt.plot(x_sel_final[i_sel][-1],y_sel_final[i_sel][-1], color = 'k', marker=cell_marker_i,
     linewidth=20, markersize=40*multiplying_factor) # end point number
    
plt.imshow(imgGray, cmap='gray')  
plt.show()


# Read in cell centers coordinates and store in cell_centers[cell][time][coordinate]
cell_centers = []
for i in range(1,201):
    with open("Mosaic_corrected/cell_" + str(i) + ".txt") as file_name:
        cell_centers.append(np.loadtxt(file_name, delimiter=","))

for cell in range(200):
    for time in range(1000):
        for coordinate in range(2):
            cell_centers[cell][time][coordinate] = 2*cell_centers[cell][time][coordinate]

# Rotate cell centers and store in rotated_cell_centers (rotated around ellipse center)
rotated_cell_centers = np.empty((200, 1000, 4))
theta_corr = np.empty(len(theta))

for i in range(1,201): # cells
  
    for j in range(len(cell_centers[0])): # time / MCS
         
        # correct theta if it shows a jump (= orientation error)  
        if theta[j]*180/np.pi >64 : # gastruloid angle
            theta_corr[j] = theta[j] - 1.57079633
        elif theta[j]*180/np.pi <0 : 
            theta_corr[j] = theta[j] + 1.57079633
        else:
            theta_corr[j] = theta[j]

        # Calculate rotated coordinates, rotate around center of ellips
        # xc/yc = center coordinates of ellipse, O_x/O_y = coordinates of cell
        O_x = np.cos(theta_corr[j]) * (cell_centers[i-1][j][0]-xc[j]) + np.sin(theta_corr[j]) * (cell_centers[i-1][j][1]-yc[j])
        rotated_cell_centers[i-1,j,0] = O_x + xc[j] # rotated x-coordinate of cell
        rotated_cell_centers[i-1,j,2] = O_x # x component of vector from center of ellipse to position of cell
        O_y = np.cos(theta_corr[j]) * (cell_centers[i-1][j][1]-yc[j]) - np.sin(theta_corr[j]) * (cell_centers[i-1][j][0]-xc[j])
        rotated_cell_centers[i-1,j,1] = O_y + yc[j] # rotated y-coordinate of cell
        rotated_cell_centers[i-1,j,3] = O_y # y component of vector from center of ellipse to position of cell


# Investigate theta_corr
%matplotlib inline
print('Theta at t_start: ', theta[t_start]*180/np.pi, 'and theta at t_end: ',theta[999]*180/np.pi)
print('Corrected theta at t_start: ', theta_corr[t_start]*180/np.pi, 'and corrected theta at t_end: ',theta_corr[999]*180/np.pi)
theta_x = range(0,1000,1)
# Plot theta
plt.scatter(theta_x,theta*180/np.pi, s = 1, color = 'red')
# Plot theta_corr
plt.scatter(theta_x,theta_corr*180/np.pi, s = 1, color = 'black')

# Plot t_start
plt.scatter(199,theta_corr[199]*180/np.pi, s = 30, color = 'limegreen')
plt.xlim(-10, 1010)

Theta at t_start:  127.77361178960888 and theta at t_end:  53.37818613761889
Corrected theta at t_start:  37.773611605969975 and corrected theta at t_end:  53.37818613761889

(-10.0, 1010.0)


# CORRECTION FOR DRIFT AND FOR ROTATION
# Rotate rotated cell centers back to orientation gastruloid at final MCS, plot trajectories at final MCS and correct for drift
# Plot figure: each trajectory a different color (max. 8 colors), color in the color-blind-friendly palette

%matplotlib inline

# Select image (time point) and adjust size (optional)
imgplot= plt.imread('Mosaic_corrected/image0099900.png') # final time point
plt.rcParams['figure.figsize'] = [200, 100]
# Make image gray (so no dark blue cells)
R, G, B = imgplot[:,:,0], imgplot[:,:,1], imgplot[:,:,2]
imgGray = 0.1 * R + 0.1 * G + 0.7*B #0.2989 * R + 0.5870 * G + 0.1140 * B

# Calculate drift of gastruloid for each MCS w.r.t. final MCS using the difference in position of ellipse centers
x_corr_drift_all = np.empty((1000))
y_corr_drift_all = np.empty((1000))
# Rotate all rotated cell centers back to gastruloid's orientation at the final MCS (999)
rotated_cell_centers_rot = np.empty((200, 1000, 4))

cmap = [#'#000000', #black
        '#e69d00', #orange
        '#56b4e9', #sky blue
        '#009e73', #bluish green
        '#f0e442', #yellow
        '#0072b2', #blue
        '#d55e00', #vermillion (reddish)
        '#7f7f7f', #grey
        '#cc79a7']#, #reddish purple
        #'#ffffff', #white
        #'#7f7f7f'] #grey
cell_marker = ['$1$','$2$','$3$','$4$','$5$','$6$','$7$','$8$','$9$','$10$','$11$','$12$','$13$','$14$','$15$','$16$','$17$','$18$','$19$','$20$']
markers_numbered_selected = [11,10,14,17,16,18,3,5,6,7,8,9]

for i in range(0,200): # cells
    
    for j in range(0,1000): # MCSs
        x_corr_drift_all[j] = xc[-1]-xc[j]
        y_corr_drift_all[j] = yc[-1]-yc[j]
        
        O_x_rot = np.cos(2*np.pi - theta_corr[-1]) * (rotated_cell_centers[i][j][0]-xc[j]) + np.sin(2*np.pi - theta_corr[-1]) * (rotated_cell_centers[i][j][1]-yc[j])
        rotated_cell_centers_rot[i,j,0] = O_x_rot + xc[j] # rotated x-coordinate of cell
        rotated_cell_centers_rot[i,j,2] = O_x_rot # x component of vector from center of ellipse to position of cell
        
        O_y_rot = np.cos(2*np.pi - theta_corr[-1]) * (rotated_cell_centers[i][j][1]-yc[j]) - np.sin(2*np.pi - theta_corr[-1]) * (rotated_cell_centers[i][j][0]-xc[j])
        rotated_cell_centers_rot[i,j,1] = O_y_rot + yc[j] # rotated y-coordinate of cell
        rotated_cell_centers_rot[i,j,3] = O_y_rot # y component of vector from center of ellipse to position of cell

for i in range(len(markers_numbered_selected)): # loop over cells
    i_sel_rot = (markers_numbered_selected[i])*10-1
    i_sel = (markers_numbered_selected[i])-1
    cell_marker_i = cell_marker[i_sel]
    
    if i == 0:
        color_i = cmap[2]
    if i == 2:
        color_i = cmap[7]
    if i > 0 and i < 8 and i != 2:
        color_i = cmap[i]
    if i == 8:
        color_i = cmap[0]
    if i == 9:
        color_i = cmap[2]
    if i == 10:
        color_i = cmap[1]
    if i == 11:
        color_i = cmap[3]
    
    plt.plot(rotated_cell_centers_rot[i_sel_rot,selected_points,0] + x_corr_drift_all[selected_points], rotated_cell_centers_rot[i_sel_rot,selected_points,1] + y_corr_drift_all[selected_points], color = 'k',
             linewidth=25, markersize=25) # line connecting start and end point black
    plt.plot(rotated_cell_centers_rot[i_sel_rot,selected_points,0] + x_corr_drift_all[selected_points], rotated_cell_centers_rot[i_sel_rot,selected_points,1] + y_corr_drift_all[selected_points], color = color_i,
             linewidth=20, markersize=20) # line connecting start and end point, linewidth=25, markersize=25

for i in range(len(markers_numbered_selected)):
    
    if i <= 5:
        multiplying_factor = 1.5 # use a larger icon for cells that start to count from 10
    else:
        multiplying_factor = 1
    
    i_sel_rot = (markers_numbered_selected[i])*10-1
    i_sel = markers_numbered_selected[i]-1 #i
    cell_marker_i = cell_marker[i_sel]
    if i == 0:
        color_i = cmap[2]
    if i == 2:
        color_i = cmap[7]
    if i > 0 and i < 8 and i != 2:
        color_i = cmap[i]
    if i == 8:
        color_i = cmap[0]
    if i == 9:
        color_i = cmap[2]
    if i == 10:
        color_i = cmap[1]
    if i == 11:
        color_i = cmap[3]
    
    # Calculate drift for starting point w.r.t. the final point and add that to starting positions
    x_corr_drift_f = xc[999]-xc[t_start]
    y_corr_drift_f = yc[999]-yc[t_start]
    # plot start point
    plt.plot(rotated_cell_centers_rot[i_sel_rot,selected_points[0],0] + x_corr_drift_all[selected_points[0]], rotated_cell_centers_rot[i_sel_rot,selected_points[0],1] + y_corr_drift_all[selected_points[0]],
     color = 'k', marker='o',linewidth=2, markersize=130) # start point black edge
    plt.plot(rotated_cell_centers_rot[i_sel_rot,selected_points[0],0] + x_corr_drift_all[selected_points[0]], rotated_cell_centers_rot[i_sel_rot,selected_points[0],1] + y_corr_drift_all[selected_points[0]],
     color = color_i, marker='o', linewidth=2, markersize=120) # start point colored,  use alpha = .5 for transparency
    plt.plot(rotated_cell_centers_rot[i_sel_rot,selected_points[0],0] + x_corr_drift_all[selected_points[0]], rotated_cell_centers_rot[i_sel_rot,selected_points[0],1] + y_corr_drift_all[selected_points[0]],
     color = 'w', marker='o', linewidth=2, markersize=80)#70) # start point white
    plt.plot(rotated_cell_centers_rot[i_sel_rot,selected_points[0],0] + x_corr_drift_all[selected_points[0]], rotated_cell_centers_rot[i_sel_rot,selected_points[0],1] + y_corr_drift_all[selected_points[0]],
     color = 'k', marker=cell_marker_i, linewidth=20, markersize=40*multiplying_factor) # start point number
    
    # plot end point
    plt.plot(rotated_cell_centers_rot[i_sel_rot,selected_points[-1],0] + x_corr_drift_all[selected_points[-1]], rotated_cell_centers_rot[i_sel_rot,selected_points[-1],1] + y_corr_drift_all[selected_points[-1]],
     color = 'k', marker='D', linewidth=20, markersize=130) # end point black edge
    plt.plot(rotated_cell_centers_rot[i_sel_rot,selected_points[-1],0] + x_corr_drift_all[selected_points[-1]], rotated_cell_centers_rot[i_sel_rot,selected_points[-1],1] + y_corr_drift_all[selected_points[-1]],
     color = color_i, marker='D', linewidth=20, markersize=120) # end point colored
    plt.plot(rotated_cell_centers_rot[i_sel_rot,selected_points[-1],0] + x_corr_drift_all[selected_points[-1]], rotated_cell_centers_rot[i_sel_rot,selected_points[-1],1] + y_corr_drift_all[selected_points[-1]],
     color = 'w', marker='D', linewidth=20, markersize=80)#70) # end point white
    plt.plot(rotated_cell_centers_rot[i_sel_rot,selected_points[-1],0] + x_corr_drift_all[selected_points[-1]], rotated_cell_centers_rot[i_sel_rot,selected_points[-1],1] + y_corr_drift_all[selected_points[-1]],
     color = 'k', marker=cell_marker_i,linewidth=20, markersize=40*multiplying_factor) # end point number

plt.imshow(imgGray, cmap='gray')
#plt.savefig("Mosaic_corrected/Figure_in_silico_mosaic_corrected_cdrift_crotation_12trajectories.png", bbox_inches='tight',pad_inches = 0)

plt.show()


# Define number of cells and the selected time points you want to plot the trajectory for
No_cells = 200 # number of cells that are being tracked

# Define list of arrays for x- and y-coordinates for each tracked cell
path_x_sim_all = [None] * (No_cells)
path_y_sim_all = [None] * (No_cells)
path_x_sim_all_final = [None] * No_cells
path_y_sim_all_final = [None] * No_cells


# Store x- and y-coordinates for each tracked cell in lists 'path_x_sim_all_final' and 'path_y_sim_all_final'
for i in range(1,int(No_cells)+1,1):
    filepath = ""
    file_dir = ['Mosaic_corrected/cell_','.txt']
    filepath = file_dir[0]+str(i)+file_dir[1]
    with open(filepath) as f:
            for line in f:
                currentline = line.split(",")
                x_line = float(currentline[0])
                y_line = float(currentline[1])
                ii = int(i)-1
                path_x_sim_all[ii] = np.append(path_x_sim_all[ii],x_line)
                path_y_sim_all[ii] = np.append(path_y_sim_all[ii],y_line)
    path_x_sim_all_final[ii] = path_x_sim_all[ii][1:(len(path_x_sim_all[ii])+1)]
    path_y_sim_all_final[ii] = path_y_sim_all[ii][1:(len(path_y_sim_all[ii])+1)]


# Store x- and y-coordinates of each tracked cell for only the selected time points in lists 'x_sel_final' and 'y_sel_final'
x_sel_all = [None] * No_cells
y_sel_all = [None] * No_cells
x_sel_all_final = [None] * No_cells
y_sel_all_final = [None] * No_cells

for j in range(No_cells):
    for k in range(len(selected_points)):
        x_sel_all[j] = np.append(x_sel_all[j],2*path_x_sim_all_final[j][selected_points[k]]) # this coding can probably be improved 
        y_sel_all[j] = np.append(y_sel_all[j],2*path_y_sim_all_final[j][selected_points[k]])
    x_sel_all_final[j] = x_sel_all[j][1:(len(selected_points)+1)]
    y_sel_all_final[j] = y_sel_all[j][1:(len(selected_points)+1)]


# Color indicates orientation of direction vector
import matplotlib.colors as colors
import matplotlib.cm as cmx
import matplotlib as mpl
name = 'twilight'# or 'hsv'
cmap=mpl.colormaps[name]
#cmap = plt.cm.jet
cNorm  = colors.Normalize(vmin=-np.pi, vmax=np.pi)
scalarMap = cmx.ScalarMappable(norm=cNorm,cmap=cmap)

%matplotlib inline
imgplot_999= plt.imread('Mosaic_corrected/image0099900.png') 
plt.rcParams['figure.figsize'] = [200, 100]

# Make image gray (so no dark blue cells)
R_999, G_999, B_999 = imgplot_999[:,:,0], imgplot_999[:,:,1], imgplot_999[:,:,2]
imgGray_999 = 0.1 * R_999 + 0.1 * G_999 + 0.7*B_999

for i in range(No_cells):
    
    # Rotated and rotated back, vector downscaled by factor 20
    x_dx_rot_rot_drift = ((rotated_cell_centers_rot[i,t_end,0] - rotated_cell_centers_rot[i,t_start,0]) - (xc[t_end] - xc[t_start]))/20
    y_dy_rot_rot_drift = ((rotated_cell_centers_rot[i,t_end,1] - rotated_cell_centers_rot[i,t_start,1]) - (yc[t_end] - yc[t_start]))/20
    
    x_arrow_base_rot_rot = rotated_cell_centers_rot[i,t_end,0] - x_dx_rot_rot_drift
    y_arrow_base_rot_rot = rotated_cell_centers_rot[i,t_end,1] - y_dy_rot_rot_drift
    
    angle_cell = np.arctan2((y_dy_rot_rot_drift),(x_dx_rot_rot_drift))
    colorVal = scalarMap.to_rgba(angle_cell)
    plt.arrow(x_arrow_base_rot_rot, y_arrow_base_rot_rot, x_dx_rot_rot_drift, y_dy_rot_rot_drift, width = 1, head_width=6, head_length=6, color=colorVal)
    plt.arrow(x_arrow_base_rot_rot+300, y_arrow_base_rot_rot-290, x_dx_rot_rot_drift, y_dy_rot_rot_drift, width = 1, head_width=6, head_length=6, color=colorVal)
    
# Save figure
plt.imshow(imgGray_999, cmap='gray', alpha = .5)
#plt.savefig("Mosaic_corrected/Figure_in_silico_mosaic_corrected_t999_all_vectors_corr_drift_rot_DIV20_angle_twilight.png", bbox_inches='tight',pad_inches = 0)
plt.show()

# https://stackoverflow.com/questions/18748328/plotting-arrows-with-different-color-in-matplotlib
# use https://stackoverflow.com/questions/19576495/color-matplotlib-quiver-field-according-to-magnitude-and-direction


# Plot color wheel
%matplotlib inline
from matplotlib import cm
import matplotlib as mpl

fig = plt.figure()

display_axes = fig.add_axes([0.1,0.1,0.8,0.8], projection='polar')
norm = mpl.colors.Normalize(0,2*np.pi)
quant_steps = 2056
cb = mpl.colorbar.ColorbarBase(display_axes, cmap=cm.get_cmap('twilight_shifted',quant_steps),norm=norm, orientation='horizontal')

# aesthetics - get rid of border and axis labels                                   
cb.outline.set_visible(False)                                 
display_axes.set_axis_off()
plt.savefig("Mosaic_corrected/color_wheel_twilight.png", bbox_inches='tight',pad_inches = 0)
plt.show()

/var/folders/9w/72zthbm57f33y8wq5zqg5kjw0000gn/T/ipykernel_1167/1637029653.py:12: MatplotlibDeprecationWarning: Auto-removal of grids by pcolor() and pcolormesh() is deprecated since 3.5 and will be removed two minor releases later; please call grid(False) first.
  cb = mpl.colorbar.ColorbarBase(display_axes, cmap=cm.get_cmap('twilight_shifted',quant_steps),norm=norm, orientation='horizontal')


# Determine the length of the simulated gastruloid at t_end = 2 * long arm of ellipse fit. The long arm should be stored in the
# parameter 'b', and the short arm in parameter 'a', but the fitting tool sometimes stored them incorrectly.
# Therefore, check the value of a and b first to determine which parameter corresponds to the long or short arm.
# In our in silico example, the long arm of the fitted ellipses is on the x-axis. The short arm on the y-axis. 
# If a and b are incorrectly switched, the values stored in the y coordinates correspond to the long arm, and x coordinates correspond to the short arm.

if a[t_end] > b[t_end]:
    print('a > b, so x-axis is long axis')
    length_gastruloid = 2*a[t_end] # used to normalize the position of the cells along the long axis
    
    end_pos_on_long_axis = rotated_cell_centers[:,t_end,0] # rotated x-coordinate at t = t_end
    end_pos_on_long_axis_norm = (end_pos_on_long_axis - xc[t_end])/(length_gastruloid) # (rotated x-coordinate at t_end - x_center_of_ellipse at t_end) divided by length gastruloid
    
    # (rotated x-coordinate at t_end - rotated x-coordinate at t_start) - difference between center of ellipses to correct for drift
    moved_end_norm = ((end_pos_on_long_axis - rotated_cell_centers[:,t_start,0]) - (xc[t_end] - xc[t_start]) ) / length_gastruloid # divided by length of gastruloid to normalize
    
else:
    print('a < b, so y-axis is long axis')
    length_gastruloid = 2*b[t_end]
    
    # Rotated cell centers around ellipse center
    end_pos_on_long_axis = rotated_cell_centers[:,t_end,1] # rotated y-coordinate at t = t_end
    end_pos_on_long_axis_norm = (end_pos_on_long_axis - yc[t_end])/(length_gastruloid) # (rotated y-coordinate at t_end - y_center_of_ellipse at t_end) divided by length gastruloid
    
    # rotated y-coordinate at t_end - rotated y-coordinate at t_start
    moved_end_norm = ((end_pos_on_long_axis - rotated_cell_centers[:,t_start,1]) - (yc[t_end] - yc[t_start]) ) / length_gastruloid # to normalize, divide by length of gastruloid
    
# In vitro time lapse data
# See cells under "In vitro"

# -------- FIGURE ----------
# Figure parameters    
%matplotlib inline
font = {'size'   : 12}
matplotlib.rc('font', **font)
matplotlib.rcParams['font.sans-serif'] = "Arial"
matplotlib.rcParams['font.family'] = "sans-serif" 

plt.gcf().clear()
fig = plt.figure(1)
ax = fig.add_subplot(111)

plt.rcParams['figure.figsize'] = [10, 10]
plt.rcParams['xtick.labelsize'] = 'x-large'
plt.rcParams['ytick.labelsize'] = 'x-large'  
plt.xlabel('End position projected on long axis', fontsize = 16)
plt.ylabel('Movement along long axis', fontsize = 16)
plt.xlim(-0.55, 0.55)
plt.ylim(-0.51, 0.51)
plt.xticks([-0.5, -0.25, 0, 0.25, 0.5],['-L/2','-L/4','0','L/4','L/2'])
plt.yticks([-0.5, -0.25, 0, 0.25, 0.5],['-L/2','-L/4','0','L/4','L/2'])

# In vitro time lapse data
ax.scatter(end_pos_on_long_axis_iv_norm, moved_end_iv_norm, color = '#cc79a7', marker = 'o', label = 'in vitro',s = 50)

# Plot rotated cell centers around ellipse center
ax.scatter(end_pos_on_long_axis_norm, moved_end_norm, s = 5, color = 'k', label = 'in silico')
handles, labels = ax.get_legend_handles_labels()
lgd = ax.legend(handles, labels, loc='upper left', bbox_to_anchor=(1,1))
fig.set_size_inches((110/5)/2.54, (76/5)/2.54)

# Save figure
#plt.savefig("Mosaic_corrected/Figure_corr_movement_along_long_axis_legend.svg", bbox_extra_artists=([lgd]), bbox_inches='tight')

plt.show()

a > b, so x-axis is long axis


# Plot movement along short axis

if a[t_end] > b[t_end]: # if a is long axis, rotated y-axis is short axis
    print('a > b, so y-axis is short axis')
    length_gastruloid_short = 2*b[t_end] # short axis
    
    # Rotated cell centers around ellipse center 
    end_pos_on_short_axis = rotated_cell_centers[:,t_end,1] # rotated y-coordinate at t = t_end
    end_pos_on_short_axis_norm = (end_pos_on_short_axis - yc[t_end])/(length_gastruloid_short) # (rotated y-coordinate at t_end - y_center_of_ellipse at t_end) divided by (short) length gastruloid
    
    # (rotated y-coordinate at t_end - rotated y-coordinate at t_start) - difference between center of ellipses to correct for drift
    moved_end_short_norm = ((end_pos_on_short_axis - rotated_cell_centers[:,t_start,1]) - (yc[t_end] - yc[t_start]) ) / length_gastruloid_short # divided by length of gastruloid to normalize
    
else:  # if b is long axis, rotated x-axis is short axis
    print('a < b, so x-axis is short axis')
    length_gastruloid_short = 2*a[t_end]
    
    # Rotated cell centers around ellipse center
    end_pos_on_short_axis = rotated_cell_centers[:,t_end,0] # rotated x-coordinate at t = t_end
    end_pos_on_short_axis_norm = (end_pos_on_short_axis - xc[t_end])/(length_gastruloid_short) # (rotated y-coordinate at t_end - y_center_of_ellipse at t_end) divided by length gastruloid
    
    # rotated x-coordinate at t_end - rotated x-coordinate at t_start
    moved_end_short_norm = ((end_pos_on_short_axis - rotated_cell_centers[:,t_start,0]) - (xc[t_end] - xc[t_start]) ) / length_gastruloid_short # to normalize, divide by length of gastruloid
    
# In vitro time lapse data
# See cells under "In vitro"

# -------- FIGURE ----------
# Figure parameters    
%matplotlib inline
font = {'size'   : 12}
matplotlib.rc('font', **font)
matplotlib.rcParams['font.sans-serif'] = "Arial"
matplotlib.rcParams['font.family'] = "sans-serif" 

plt.gcf().clear()
fig = plt.figure(1)
ax = fig.add_subplot(111)

plt.rcParams['figure.figsize'] = [10, 10]
plt.rcParams['xtick.labelsize'] = 'x-large'
plt.rcParams['ytick.labelsize'] = 'x-large'  
plt.xlabel('End position projected on short axis', fontsize = 16)
plt.ylabel('Movement along short axis', fontsize = 16)
plt.xlim(-1, 1)
plt.ylim(-1, 1)
plt.xticks([-1, -0.5, 0, 0.5, 1],['S','-S/2','0','S/2', 'S'])
plt.yticks([-1,-0.5, 0, 0.5, 1],['-S','-S/2','0','S/2','S'])

ax.scatter(end_pos_on_short_axis_iv_norm, moved_end_iv_short_norm, color = '#cc79a7', marker = 'o', s = 50)#label = 'in vitro',s = 50)
ax.scatter(end_pos_on_short_axis_norm, moved_end_short_norm, s = 5, color = 'k')#, label = 'in silico')

fig.set_size_inches((110/5)/2.54, (76/5)/2.54)

# Save figure
#plt.savefig("Mosaic_corrected/Figure_corr_movement_along_short_axis_end_pos_short_eqaxS.svg", bbox_extra_artists=([lgd]), bbox_inches='tight')

plt.show()

a > b, so y-axis is short axis


# Plot movement along short axis

if a[t_end] > b[t_end]: # if a is long axis, rotated y-axis is short axis
    print('a > b, so y-axis is short axis')
    length_gastruloid_short = 2*b[t_end] # short axis
    
    # Rotated cell centers around ellipse center 
    end_pos_on_short_axis = rotated_cell_centers[:,t_end,1] # rotated y-coordinate at t = t_end
    end_pos_on_short_axis_norm = (end_pos_on_short_axis - yc[t_end])/(length_gastruloid_short) # (rotated y-coordinate at t_end - y_center_of_ellipse at t_end) divided by (short) length gastruloid
    
    # (rotated y-coordinate at t_end - rotated y-coordinate at t_start) - difference between center of ellipses to correct for drift
    moved_end_short_norm = ((end_pos_on_short_axis - rotated_cell_centers[:,t_start,1]) - (yc[t_end] - yc[t_start]) ) / length_gastruloid_short # divided by length of gastruloid to normalize
    
else:  # if b is long axis, rotated x-axis is short axis
    print('a < b, so x-axis is short axis')
    length_gastruloid_short = 2*a[t_end]
    
    # Rotated cell centers around ellipse center
    end_pos_on_short_axis = rotated_cell_centers[:,t_end,0] # rotated x-coordinate at t = t_end
    end_pos_on_short_axis_norm = (end_pos_on_short_axis - xc[t_end])/(length_gastruloid_short) # (rotated y-coordinate at t_end - y_center_of_ellipse at t_end) divided by length gastruloid
    
    # rotated x-coordinate at t_end - rotated x-coordinate at t_start
    moved_end_short_norm = ((end_pos_on_short_axis - rotated_cell_centers[:,t_start,0]) - (xc[t_end] - xc[t_start]) ) / length_gastruloid_short # to normalize, divide by length of gastruloid
    
inw_outw = np.arange(0)

for i in range(0,len(end_pos_on_short_axis_norm)):
    #print(i)
    if end_pos_on_short_axis_norm[i] > 0 and moved_end_short_norm[i] < 0:
        inw_outw = np.append(inw_outw, 1) # inwards
    elif end_pos_on_short_axis_norm[i] > 0 and moved_end_short_norm[i] > 0 and moved_end_short_norm[i] > end_pos_on_short_axis_norm[i]:
        inw_outw = np.append(inw_outw, 2) # inwards, outwards
    elif end_pos_on_short_axis_norm[i] > 0 and moved_end_short_norm[i] > 0 and moved_end_short_norm[i] < end_pos_on_short_axis_norm[i]:
        inw_outw = np.append(inw_outw, 0) # outwards
    elif end_pos_on_short_axis_norm[i] < 0 and moved_end_short_norm[i] < 0 and moved_end_short_norm[i] < end_pos_on_short_axis_norm[i]:
        inw_outw = np.append(inw_outw, 2) # inwards, outwards
    elif end_pos_on_short_axis_norm[i] < 0 and moved_end_short_norm[i] < 0 and moved_end_short_norm[i] > end_pos_on_short_axis_norm[i]:
        inw_outw = np.append(inw_outw, 0) # outwards
    elif end_pos_on_short_axis_norm[i] < 0 and moved_end_short_norm[i] > 0:
        inw_outw = np.append(inw_outw, 1) # inwards
    else:
        inw_outw = np.append(inw_outw, 3) # something is zero

# In vitro time lapse data
# See cells under "In vitro"

# -------- FIGURE ----------
# Figure parameters    
%matplotlib inline
font = {'size'   : 12}
matplotlib.rc('font', **font)
matplotlib.rcParams['font.sans-serif'] = "Arial"
matplotlib.rcParams['font.family'] = "sans-serif" 

plt.gcf().clear()
fig = plt.figure(1)
ax = fig.add_subplot(111)

plt.rcParams['figure.figsize'] = [10, 10]
plt.rcParams['xtick.labelsize'] = 'x-large'
plt.rcParams['ytick.labelsize'] = 'x-large'  
plt.xlabel('End position projected on short axis', fontsize = 16)
plt.ylabel('Movement along short axis', fontsize = 16)
plt.xlim(-1, 1)
plt.ylim(-1, 1)
plt.xticks([-1, -0.5, 0, 0.5, 1],['S','-S/2','0','S/2', 'S'])
plt.yticks([-1,-0.5, 0, 0.5, 1],['-S','-S/2','0','S/2','S'])

# In vitro time lapse data
inw_outw_iv = np.arange(0)
for i in range(0,len(end_pos_on_short_axis_iv_norm)):
    #print(i)
    if end_pos_on_short_axis_iv_norm[i] > 0 and moved_end_iv_short_norm[i] < 0:
        inw_outw_iv = np.append(inw_outw_iv, 1) # inwards
    elif end_pos_on_short_axis_iv_norm[i] > 0 and moved_end_iv_short_norm[i] > 0 and moved_end_iv_short_norm[i] > end_pos_on_short_axis_iv_norm[i]:
        inw_outw_iv = np.append(inw_outw_iv, 2) # inwards, outwards
    elif end_pos_on_short_axis_iv_norm[i] > 0 and moved_end_iv_short_norm[i] > 0 and moved_end_iv_short_norm[i] < end_pos_on_short_axis_iv_norm[i]:
        inw_outw_iv = np.append(inw_outw_iv, 0) # outwards
    elif end_pos_on_short_axis_iv_norm[i] < 0 and moved_end_iv_short_norm[i] < 0 and moved_end_iv_short_norm[i] < end_pos_on_short_axis_iv_norm[i]:
        inw_outw_iv = np.append(inw_outw_iv, 2) # inwards, outwards
    elif end_pos_on_short_axis_iv_norm[i] < 0 and moved_end_iv_short_norm[i] < 0 and moved_end_iv_short_norm[i] > end_pos_on_short_axis_iv_norm[i]:
        inw_outw_iv = np.append(inw_outw_iv, 0) # outwards
    elif end_pos_on_short_axis_iv_norm[i] < 0 and moved_end_iv_short_norm[i] > 0:
        inw_outw_iv = np.append(inw_outw_iv, 1) # inwards
    else:
        inw_outw_iv = np.append(inw_outw_iv, 3) # something is zero


color_i = ['blue', #outwards
           'orange', #inwards
           'skyblue', #inwards outwards
           'yellow'] # not sure]

for i in range(0,len(end_pos_on_long_axis_iv_norm)):
    ax.scatter(end_pos_on_short_axis_iv_norm[i], moved_end_iv_short_norm[i], color = color_i[inw_outw_iv[i]], marker = 'o',s = 50)

for i in range(0,len(end_pos_on_short_axis_norm)):
    ax.scatter(end_pos_on_short_axis_norm[i], moved_end_short_norm[i], s = 5, color = color_i[inw_outw[i]])#, label = 'in silico')

fig.set_size_inches((110/5)/2.54, (76/5)/2.54)

# Save figure
#plt.savefig("Mosaic_corrected/Figure_corr_movement_along_short_axis_end_pos_short_colored_eqaxS.svg", bbox_extra_artists=([lgd]), bbox_inches='tight')

plt.show()

a > b, so y-axis is short axis


# Plot movement along short axis vs. end position projected on long axis

if a[t_end] > b[t_end]: # if a is long axis, rotated y-axis is short axis
    print('a > b, so y-axis is short axis')
    length_gastruloid_short = 2*b[t_end] # short axis
    length_gastruloid_long = 2*a[t_end]
    # Rotated cell centers around ellipse center 
    end_pos_on_short_axis = rotated_cell_centers[:,t_end,1] # rotated y-coordinate at t = t_end
    end_pos_on_long_axis = rotated_cell_centers[:,t_end,0] # rotated x-coordinate at t = t_end
    end_pos_on_short_axis_norm = (end_pos_on_short_axis - yc[t_end])/(length_gastruloid_short) # (rotated y-coordinate at t_end - y_center_of_ellipse at t_end) divided by (short) length gastruloid
    end_pos_on_long_axis_norm = (end_pos_on_long_axis - xc[t_end])/(length_gastruloid_long) # (rotated x-coordinate at t_end - y_center_of_ellipse at t_end) divided by (short) length gastruloid
    
    # (rotated y-coordinate at t_end - rotated y-coordinate at t_start) - difference between center of ellipses to correct for drift
    moved_end_short_norm = ((end_pos_on_short_axis - rotated_cell_centers[:,t_start,1]) - (yc[t_end] - yc[t_start]) ) / length_gastruloid_short # divided by length of gastruloid to normalize
    moved_end_long_norm = ((end_pos_on_long_axis - rotated_cell_centers[:,t_start,0]) - (xc[t_end] - xc[t_start]) ) / length_gastruloid_long # divided by length of gastruloid to normalize

else:  # if b is long axis, rotated x-axis is short axis
    print('a < b, so x-axis is short axis, y-axis is long axis')
    length_gastruloid_short = 2*a[t_end]
    length_gastruloid_long = 2*b[t_end]
    
    # Rotated cell centers around ellipse center
    end_pos_on_short_axis = rotated_cell_centers[:,t_end,0] # rotated x-coordinate at t = t_end
    end_pos_on_long_axis = rotated_cell_centers[:,t_end,1] # rotated y-coordinate at t = t_end
    end_pos_on_short_axis_norm = (end_pos_on_short_axis - xc[t_end])/(length_gastruloid_short) # (rotated y-coordinate at t_end - y_center_of_ellipse at t_end) divided by length gastruloid
    end_pos_on_long_axis_norm = (end_pos_on_long_axis - yc[t_end])/(length_gastruloid_long) # (rotated x-coordinate at t_end - y_center_of_ellipse at t_end) divided by (short) length gastruloid
    
    # rotated x-coordinate at t_end - rotated x-coordinate at t_start
    moved_end_short_norm = ((end_pos_on_short_axis - rotated_cell_centers[:,t_start,0]) - (xc[t_end] - xc[t_start]) ) / length_gastruloid_short # to normalize, divide by length of gastruloid
    moved_end_long_norm = ((end_pos_on_long_axis - rotated_cell_centers[:,t_start,1]) - (yc[t_end] - yc[t_start]) ) / length_gastruloid_long # divided by length of gastruloid to normalize

inw_outw = np.arange(0)

for i in range(0,len(end_pos_on_short_axis_norm)):
    #print(i)
    if end_pos_on_short_axis_norm[i] > 0 and moved_end_short_norm[i] < 0:
        inw_outw = np.append(inw_outw, 1) # inwards
    elif end_pos_on_short_axis_norm[i] > 0 and moved_end_short_norm[i] > 0 and moved_end_short_norm[i] > end_pos_on_short_axis_norm[i]:
        inw_outw = np.append(inw_outw, 2) # inwards, outwards
    elif end_pos_on_short_axis_norm[i] > 0 and moved_end_short_norm[i] > 0 and moved_end_short_norm[i] < end_pos_on_short_axis_norm[i]:
        inw_outw = np.append(inw_outw, 0) # outwards
    elif end_pos_on_short_axis_norm[i] < 0 and moved_end_short_norm[i] < 0 and moved_end_short_norm[i] < end_pos_on_short_axis_norm[i]:
        inw_outw = np.append(inw_outw, 2) # inwards, outwards
    elif end_pos_on_short_axis_norm[i] < 0 and moved_end_short_norm[i] < 0 and moved_end_short_norm[i] > end_pos_on_short_axis_norm[i]:
        inw_outw = np.append(inw_outw, 0) # outwards
    elif end_pos_on_short_axis_norm[i] < 0 and moved_end_short_norm[i] > 0:
        inw_outw = np.append(inw_outw, 1) # inwards
    else:
        inw_outw = np.append(inw_outw, 3) # something is zero
    
# In vitro time lapse data
# See cells under "In vitro"
inw_outw_iv = np.arange(0)
for i in range(0,len(end_pos_on_short_axis_iv_norm)):
    #print(i)
    if end_pos_on_short_axis_iv_norm[i] > 0 and moved_end_iv_short_norm[i] < 0:
        inw_outw_iv = np.append(inw_outw_iv, 1) # inwards
    elif end_pos_on_short_axis_iv_norm[i] > 0 and moved_end_iv_short_norm[i] > 0 and moved_end_iv_short_norm[i] > end_pos_on_short_axis_iv_norm[i]:
        inw_outw_iv = np.append(inw_outw_iv, 2) # inwards, outwards
    elif end_pos_on_short_axis_iv_norm[i] > 0 and moved_end_iv_short_norm[i] > 0 and moved_end_iv_short_norm[i] < end_pos_on_short_axis_iv_norm[i]:
        inw_outw_iv = np.append(inw_outw_iv, 0) # outwards
    elif end_pos_on_short_axis_iv_norm[i] < 0 and moved_end_iv_short_norm[i] < 0 and moved_end_iv_short_norm[i] < end_pos_on_short_axis_iv_norm[i]:
        inw_outw_iv = np.append(inw_outw_iv, 2) # inwards, outwards
    elif end_pos_on_short_axis_iv_norm[i] < 0 and moved_end_iv_short_norm[i] < 0 and moved_end_iv_short_norm[i] > end_pos_on_short_axis_iv_norm[i]:
        inw_outw_iv = np.append(inw_outw_iv, 0) # outwards
    elif end_pos_on_short_axis_iv_norm[i] < 0 and moved_end_iv_short_norm[i] > 0:
        inw_outw_iv = np.append(inw_outw_iv, 1) # inwards
    else:
        inw_outw_iv = np.append(inw_outw_iv, 3) # something is zero

# -------- FIGURE ----------
# Figure parameters    
%matplotlib inline
font = {'size'   : 12}
matplotlib.rc('font', **font)
matplotlib.rcParams['font.sans-serif'] = "Arial"
matplotlib.rcParams['font.family'] = "sans-serif" 

plt.gcf().clear()
fig = plt.figure(1)
ax = fig.add_subplot(111)

plt.rcParams['figure.figsize'] = [10, 10]
plt.rcParams['xtick.labelsize'] = 'x-large'
plt.rcParams['ytick.labelsize'] = 'x-large'  
plt.xlabel('End position projected on long axis', fontsize = 16)
plt.ylabel('Movement along short axis', fontsize = 16)
plt.xlim(-0.55, 0.55)
plt.ylim(-1, 1)#-0.38, 0.38)
plt.xticks([-0.5, -0.25, 0, 0.25, 0.5],['-L/2','-L/4','0','L/4','L/2'])
plt.yticks([-1,-0.5, 0, 0.5,1],['-S','-S/2','0','S/2','S'])

color_i = ['blue', #outwards
           'orange', #inwards
           'skyblue', #inwards outwards
           'yellow'] # not sure]

# In vitro time lapse data
for i in range(0,len(end_pos_on_long_axis_iv_norm)):
    ax.scatter(end_pos_on_long_axis_iv_norm[i], moved_end_iv_short_norm[i], color = color_i[inw_outw_iv[i]], marker = 'o',s = 50)

# Plot rotated cell centers around ellipse center
for i in range(0,len(end_pos_on_short_axis_norm)):
    ax.scatter(end_pos_on_long_axis_norm[i], moved_end_short_norm[i], s = 5, color = color_i[inw_outw[i]])#, label = 'in silico')

fig.set_size_inches((110/5)/2.54, (76/5)/2.54)

# Save figure
#plt.savefig("Mosaic_corrected/Figure_corr_movement_along_short_axis_end_pos_long_colored_more_spaceS.svg", bbox_extra_artists=([lgd]), bbox_inches='tight')

plt.show()

a > b, so y-axis is short axis


# First, change image type from 16 bit (raw tiff file) or 32 bit (prediction) to 8 bit, e.g. in FIJI.

# Read RGB images
%matplotlib inline

# Read raw and predicted image
#image_raw = cv2.imread('Mosaic_TL_mCherry/N2V_comp/T1-z11-raw-crop_8bit.tif',0)
#image_pred = cv2.imread('Mosaic_TL_mCherry/N2V_comp/prediction-t1-200epochs-z11_8bit_crop.tif',0)

#image_raw = cv2.imread('Mosaic_TL_mCherry/N2V_comp/t5-z10-crop_8bit.tif',0)
#image_pred = cv2.imread('Mosaic_TL_mCherry/N2V_comp/prediction-t5_t7-model-200epochs-z10-crop_8bit.tif',0)

#image_raw = cv2.imread('Mosaic_TL_mCherry/N2V_comp/t23-z6-crop_8bit.tif',0)
#image_pred = cv2.imread('Mosaic_TL_mCherry/N2V_comp/prediction-t23_t31-model-200epochs_z6_crop_8bit.tif',0)

image_raw = cv2.imread('Mosaic_TL_mCherry/N2V_comp/MAX_T31_crop-raw_8bit.tif',0)
image_pred = cv2.imread('Mosaic_TL_mCherry/N2V_comp/MAX_prediction-t31-200epochs-crop_8bit.tif',0)

plt.imshow(image_raw,cmap='gray')
plt.show()
plt.imshow(image_pred,cmap='gray')
plt.show()


# create the histogram
counts_raw, bins_raw = np.histogram(image_raw, bins=256, range=(0, 255))
counts_pred, bins_pred = np.histogram(image_pred, bins=256, range=(0, 255))

%matplotlib inline
plt.bar(bins_raw[:-1],counts_raw,width=1)
plt.bar(bins_pred[:-1],counts_pred,width=1)

<BarContainer object of 256 artists>


# flatten the 2D image arrays
image_raw_flat = image_raw.flatten()
image_pred_flat = image_pred.flatten()

# Calculate 5th percentile of raw data, such that that value can be used as the minimum displayed value
image_raw_sort = np.sort(image_raw_flat)
print(image_raw_sort)
q5_raw=int(np.round(0.05*(len(image_raw_flat))))
print(q5_raw)
print('minimum displayed value for raw 8-bit image is: ', image_raw_sort[q5_raw])

# Calculate 5th percentile of prediction data, such that that value can be used as the minimum displayed value
image_pred_sort = np.sort(image_pred_flat)
print(image_pred_sort)
q5_pred=int(np.round(0.05*(len(image_pred_flat))))
print(q5_pred)
print('minimum displayed value for pred. 8-bit image is: ', image_pred_sort[q5_pred])

[  2   3   3 ... 255 255 255]
22714
minimum displayed value for raw 8-bit image is:  16
[  1   1   1 ... 255 255 255]
22714
minimum displayed value for pred. 8-bit image is:  4

A combination of convergent extension and differential adhesion explains the shapes of elongating gastruloids¶

Import the required libraries¶

LOCO-EFA of gastruloid shapes¶

Calculate the LOCO-EFA values of gastruloids shapes¶

0. Create the environment from the environment.yml file¶

1. Use our segmentation tool to create outlines from fluorescently labelled gastruloids¶

2. Store each outline into a separate csv file¶

3. Calculate for each outline file the LOCO-EFA values¶

Define which gastruloids are unique and plot their L2/L1 vs L3/L1 values (96h, 72h)¶

Place the numbers belonging to unique gastruloids in one sequence.¶

Read in the txt file with L values for each gastruloid and divide L3 and L2 by L1 (normalization)¶

Plot the figures¶

In vitro time lapse analysis¶

0a. Create binary image¶

0b. Extract outlines from binary images¶

1. Time lapse: load outlines, fit ellipse and store ellipse parameters¶

2. Cell tracking: convert ROIs to arrays and store in list¶

3. Analysis: rotate the positions of the cells to calculate their position w.r.t. the long axis¶

4. Plot: starting points of 14 cells on t = 1 image¶

5a. Plot: trajectories of 14 cells on t = 31 image¶

5b. Plot: correct trajectory plot for drift of gastruloid only¶

5c. Plot (FIG 5A): correct trajectory plot for drift and rotation of gastruloid¶

6a. Shape analysis (FIG 5B): plot cell positions on t1 + yellow connecting lines¶

6b. Shape analysis (FIG 5B): plot cell positions on t31 + yellow connecting lines¶

7a. Analysis/plot: calculate the distance of cells moved along the long axis and plot for final position on long axis¶

7b. Analysis/plot: calculate the distance moved along short axis and plot for final position on long axis¶

In silico¶

Create a binary image¶

Import images and create a GIF¶

In silico tracking - mosaic example¶

0a. Create mask and extract outline¶

0b. Load outline and fit ellipse, store ellipse parameters¶

1a. Store x- and y-coordinates of selected cells (20) for selected MCS (16) in arrays¶

2. Plot: 12 start positions on start image¶

3. Plot: and 12 end positions on end image¶

4. [Trajectories preparation] Rotate cell centers with found theta from ellipse fit, rotate around ellipse center¶

5. Plot: trajectories on end image (+ correction for drift + correction for rotation)¶

1b. OR store x- and y-coordinates of all cells (200) for selected MCS (16) in array¶

6. Plot: direction vector of each cell on final image using rotated-rotated cell centers, and correction for drift¶

7a. Plot: 'Movement along long axis' vs 'End position projected on long axis', normalized¶

7b. Plot: 'Movement along short axis' vs 'End position projected on short axis', normalized, colored¶

7c. Plot: 'Movement along short axis' vs. 'End position projected on long axis', normalized, colored¶

Noise2Void comparison¶

Calculate 5th percentile of both raw and predicted image data => use that value as the minimum displayed value¶